OSCR

Dynamic Functional Synchronization Profiles in Autism Differ by Spatial Scale and Along Hierarchical Cortical Gradients.

Code ↔ Paper

7 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 7 matches
  1. [1] § Methods and Materials › Turbulence Analysis ↔ code/turbulence.m, lines 69–167 · score 0.89 · analytic signal, Euclidean distance, Hilbert transform, kernel function, Kuramoto local, summed
  2. [2] § Methods and Materials › Statistical Analysis ↔ code/Turbulence_postproc_analysis.Rmd, lines 634–692 · score 0.73 · information cascade flow, linear models, information transfer, sex, Cohen, ADOS
  3. [3] § Methods and Materials › Turbulence Analysis ↔ code/Turbulence_postproc_analysis.Rmd, lines 634–692 · score 0.67 · standard deviation, information cascade flow, information transfer, vector, linear, row
  4. [4] § Results › Variability of Parcel-Level Functional Synchronization Maps Onto S-A and Functional Gradients in Autism ↔ code/gradient_computations.ipynb, lines 115–193 · score 0.58 · Spearman rho, turbulence maps, functional gradient, fitted, ranks, Cohen
  5. [5] § Methods and Materials › Turbulence Analysis ↔ code/turbulence.m, lines 69–167 · score 0.55 · Kuramoto local, information transfer, vector, row, dimensional, computation
  6. [6] § Results › Variability of Parcel-Level Functional Synchronization Maps Onto S-A and Functional Gradients in Autism ↔ code/gradient_computations.ipynb, lines 115–193 · score 0.55 · turbulence maps, 0.07 mm, Spearman, permutation, spin, gradients
  7. [7] § Results › Weaker, Less Spatially Correlated, and Rapidly Decaying Functional Synchronization Over Space in Autism ↔ code/Turbulence_postproc_analysis.Rmd, lines 341–442 · score 0.51 · ADOS percentage scores, information transfer, intermediate, KLOP, parcels

Paper

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

The paper is loaded when this pane is shown.

The authors' code

R Markdown · 1,459 lines · 57 KB · no license · 3 matches

  1. ---
  2. title: "Turbulence_final_revision"
  3. output: html_document
  4. date: "2026-03-21"
  5. ---
  7. Post-Processing analysis of Turbulence Project on Abide data.
  9. ```{r, message=FALSE, warning=FALSE}
  10. ### Loading libraries
  11. library("easypackages")
  12. libraries("scales","reticulate", "ppcor","gganimate", "patchwork","ggseg","ggsegSchaefer","car","gvlma","performance","effectsize", "grid","gridExtra","MatchIt","gt","lme4","lsr", "effsize", "ggsignif","tidyverse","here","sva","R.matlab","car","ggplot2", "tidyverse","ggeasy","ggpackets","emmeans", "plotly", "fmsb", "ggradar", "ggseg3d", "minpack.lm", "psych")
  14. ### Handling paths
  15. root = Sys.getenv("ROOT_PATH")
  16. #root = "/Users/bened/Library/CloudStorage/OneDrive-FondazioneIstitutoItalianoTecnologia/Desktop/University/PhD_IIT/turbulence/Turbulence_Autisms/Turbulence_Project/turbulence_final_folder" # For local processing
  17. figures_directory = paste0(root,"/figures")
  19. ### Loading functions and data
  20. #Geom_scatterbox for plotting
  21. geom_scatterbox <- ggpacket() +
  22. geom_jitter(size = 1.5, width = 0.25) +
  23. geom_boxplot(fill = NA, colour = "#000000", outlier.shape = NA, size = 1)
  25. #Loading in the phenotypical and qualitative subject and scanner data
  26. filenames = list.files(file.path(root, "data","postproc"))
  27. motionTable = read.csv(file.path(root, "data", "pheno", "qc_preproc_abideI_abideII.csv"))
  28. qcLombardo = read.csv(file.path(root, "data", "pheno", "qc_preproc_mvlombardo.csv"))
  29. mri_data = read.csv2((file.path(root, "data", "pheno", "mriparams_summary_abide.csv")))
  30. pheno = read.csv(file.path(root, "data", "pheno","pheno_abide_I_abide_II.csv"))
  32. #### Data Filtering
  33. #Filtering data for FD and Quality-checks
  34. #Merging dataframes by subid, site and dataset
  35. mergedPheno = merge.data.frame(motionTable,pheno, by = c("dataset", "site","subid"))
  36. mergedPheno = merge.data.frame(mergedPheno, qcLombardo, c("dataset", "site","subid"))
  37. mergedPheno = merge.data.frame(mergedPheno, mri_data, c("dataset", "site"))
  39. #Extracting all the IDs for the subjects that pass the quality check and have FD < 0.5
  40. filteredPheno = mergedPheno %>%
  41. dplyr::filter(pass == "YES") %>%
  42. dplyr::filter(mean_fd < 0.5)
  44. #Check visually for distributions of diagnosis over sites
  45. table(filteredPheno$site, filteredPheno$diagnosis)
  47. # merge LEUVEN, NYU, and UM labels, to even distributions in these sites
  48. filteredPheno$site_orig = filteredPheno$site
  49. filteredPheno$site[is.element(filteredPheno$site, c("NYU_1","NYU_2"))] = "NYU"
  50. filteredPheno$site[is.element(filteredPheno$site, c("LEUVEN_1","LEUVEN_2"))] = "LEUVEN"
  51. filteredPheno$site[is.element(filteredPheno$site, c("UM_1","UM_2"))] = "UM"
  52. filteredPheno$site[is.element(filteredPheno$site, c("UCLA_1","UCLA_2"))] = "UCLA"
  53. filteredPheno$site = factor(filteredPheno$site)
  55. #Check distributions again
  56. table(filteredPheno$site, filteredPheno$diagnosis)
  58. ```
  60. ```{r warning=FALSE}
  62. # remove UM, since it has an odds ratio of almost 9:1
  63. mask = filteredPheno$site != "UM"
  64. filteredPheno = filteredPheno %>% dplyr::filter(mask)
  65. filteredPheno$site = factor(filteredPheno$site)
  66. filteredPheno$mean_fd = as.numeric(filteredPheno$mean_fd)
  68. #Check, if site and diagnosis are equally distributed across the sites with a chi-square test
  69. chisq.test(filteredPheno$site, filteredPheno$diagnosis)
  70. ```
  72. p-value > 0.05, so we can go on.
  74. ```{r}
  75. # Checking distributions of framewise displacement.
  76. t.test(mean_fd ~ as.factor(diagnosis), data = filteredPheno)
  77. ```
  78. Significant (p < 0.001). Autistic subjects move more in the scanner.
  80. ```{r}
  81. #Plotting the association of FD and diagnosis
  82. t_meanfd = t.test(data = filteredPheno, mean_fd ~ diagnosis)
  83. p_meanfd = t_meanfd$p.value
  84. tab = table(filteredPheno$diagnosis)
  85. d_meanfd = round(t_meanfd$statistic / sqrt((as.double(tab[1])*as.double(tab[2]))/(as.double(tab[1])+as.double(tab[2]))),3)
  86. mean_fd_dx_plot_before = ggplot(data = filteredPheno, aes(x = diagnosis, y = mean_fd, colour = diagnosis)) +
  87. geom_scatterbox() +
  88. geom_signif(comparisons = list(c("Autism", "Control")),
  89. map_signif_level = TRUE,
  90. test = "t.test",
  91. textsize = 8,
  92. color = "black",
  93. size = 1.3) +
  94. ylab("Mean Framewise Displacement (FD) mm") +
  95. xlab("Diagnosis") + guides(colour = "none") +
  96. ylim(NA, max(filteredPheno$mean_fd) + 0.1) +
  97. theme_minimal() +
  98. theme(
  99. axis.title = element_blank(),
  100. plot.title = element_blank()
  101. ) +
  102. scale_color_manual(values = c("Control" = "#1f77b4", "Autism" = "#ff7f0e")) +
  103. ggtitle("Mean Framewise displacement of ASD vs. Control")
  105. mean_fd_dx_plot_before
  106. ```
  107. Addressing the increased movement of autism vs TD via matching the subjects from both groups with a similarity of +- 0.05mm.
  109. ```{r}
  110. #factorising the diagnosis and sex variables
  111. filteredPheno$sex = as.factor(filteredPheno$sex)
  112. filteredPheno$diagnosis = as.factor(filteredPheno$diagnosis)
  114. #matching the two groups on FD
  115. m_1 = matchit(formula = diagnosis ~ mean_fd + sex + age_years,
  116. data = filteredPheno,
  117. method = "nearest",
  118. caliper = 0.05)
  119. summary(m_1)
  120. plot(m_1,
  121. type = "density",
  122. which.xs = ~mean_fd,
  123. interactive = FALSE)
  124. ```
  127. ```{r}
  128. # Successful matching; extracting the data
  129. filteredPheno = match_data(m_1)
  130. t.test(data= filteredPheno, mean_fd ~ diagnosis)
  131. t_meanfd = t.test(data = filteredPheno, mean_fd ~ diagnosis)
  132. p_meanfd = t_meanfd$p.value
  133. d_meanfd = round(t_meanfd$statistic / sqrt((as.double(tab[1])*as.double(tab[2]))/(as.double(tab[1])+as.double(tab[2]))),3)
  134. mean_fd_dx_plot_after = ggplot(data = filteredPheno, aes(x = diagnosis, y = mean_fd, colour = diagnosis)) +
  135. geom_scatterbox() +
  136. geom_signif(comparisons = list(c("Autism", "Control")),
  137. map_signif_level = TRUE,
  138. test = "t.test",
  139. textsize = 8,
  140. color = "black",
  141. size = 1.3) +
  142. ylab("Mean Framewise Displacement (FD) mm") +
  143. xlab("Diagnosis") + guides(colour = "none") +
  144. ylim(NA, max(filteredPheno$mean_fd) + 0.1) +
  145. theme_minimal() +
  146. scale_color_manual(values = c("Control" = "#1f77b4", "Autism" = "#ff7f0e")) +
  147. ggtitle("Mean Framewise displacement of Autism vs. Control")
  149. mean_fd_dx_plot_after
  150. ```
  152. ```{r}
  153. ### Loading turbulence data
  155. #Variable pre-allocation
  156. flattenedData = c()
  157. turbData = list()
  158. allIDs = list()
  159. out_data = list()
  160. counter = 0
  161. names_turb = names(turbData)
  162. names_allIDs = as.character(allIDs$subid)
  164. ### Exclusion of Subjects based on fallible Hilbert-transormation
  165. #For 5 of the subjects, the time-series extraction did not work. The phases matrix is empty. Exclude these from the analysis.
  166. outliers_subid = c(51468, 28793, 50726, 28939, 51193)
  167. for (sub in outliers_subid){
  168. currID = as.character(sub)
  169. currName = paste0(currID,"_turbulence_measures.mat")
  170. matData = readMat(file.path(root, "data","postproc", currName))
  171. matData = matData[["output"]]
  172. matData = data.frame(matData)
  173. matData = matData$X1.1
  174. out_data[[as.character(sub)]] = matData
  175. }
  178. #Looping over the subjects of the phenotypical data to read in Turbulence measures. This step is necessary since the subids in the phenotypical are not matching with all the subids of which fMRI data was available.
  180. for (subject in 1:length(filteredPheno$subid)){
  182. #Extract the current subject-ID
  183. currID = filteredPheno$subid[subject]
  184. currName = paste0(currID,"_turbulence_measures.mat")
  186. #Skip outliers
  187. if (currID %in% outliers_subid){
  188. next}
  190. #If that ID is in the list of subjects on which the turbulence analysis was applied, go on, otherwise, exit
  191. if (currName %in% filenames){
  192. #Reading the Turbulence values of current subject
  193. matData = readMat(file.path(root, "data","postproc", currName))
  194. matData = matData[["output"]]
  195. matData = data.frame(matData)
  196. matData = matData$X1.1
  197. #Deleting the first entry of the Information.Cascade flow (because they are all NAs, due to formal, theoretical reasons)
  198. matData[[6]] = matData[[6]][-1]
  199. #For the information transfer, we still have to subtract the values calculated by the MATLAB script
  200. #from a constant value k (k=2). The reason for that is, to keep the rational consistent with the other measures:
  201. #higher values mean the information travels across longer distances (originally in the MATLAB script
  202. #higher values meant steeper slope which in turn means information travels across smaller distances)
  203. #For all that, see methods of Cruzat et al., 2022)
  204. matData[[4]] = 2 - matData[[4]]
  206. if (length(matData) == 11){
  207. #Compute the mean across time of the KLOP for each parcel
  208. matData$mean_KLOP = apply(matData$LocalKuramoto, c(1,2), FUN = mean)
  209. names(matData)[c(9,12)] = names(matData)[c(12,9)]
  210. matData[c(9, 12)] = matData[c(12, 9)]
  211. }
  213. #Delete these two, because not needed
  214. matData$Phases = NULL
  215. matData$LocalKuramoto = NULL
  217. #Some of the subjects have fallible node turbulence (only extracted for 1000 nodes). Leave them out.
  218. if (dim(matData$Turbulence.Node)[1] == 1054){
  219. counter = counter + 1
  220. allIDs[[length(allIDs)+1]] = currID #Save that ID in a list
  221. #Saving the entry in a new list
  222. turbData[[as.character(currID)]] = matData #List-Variable with all the Turbulence measures for all subjects
  223. }
  224. }
  225. }
  227. ### Performing batch correction to control for differences in scanning sites
  229. # Extracting number of parcels and lambda values from turbulence data
  230. nParcels = dim(matData$Turbulence.Node)[1]
  231. nLambda = dim(matData$Turbulence.Node)[2]
  232. dim_node = nParcels*nLambda
  233. dim_RSN = dim(matData$Turbulence.RSN)[1] * nLambda
  234. dim_output = length(matData)
  236. #Create data.frame for subject IDs
  237. allIDs = unlist(allIDs)
  238. allIDs = data.frame(allIDs)
  239. colnames(allIDs) = "subid"
  240. #Merge the analyzed subjects with the phenotypical data
  241. filteredPheno = merge.data.frame(filteredPheno,allIDs, by = "subid")
  242. #Order it in the same way as my Turbulence matrix
  243. filteredPheno = filteredPheno[match(allIDs$subid, filteredPheno$subid),]
  244. #Factorise the covariates site, sex and diagnosis and numerise mean_fd and age
  245. filteredPheno = filteredPheno %>% mutate(
  246. sex = as.factor(sex),
  247. site = as.factor(site),
  248. diagnosis = as.factor(diagnosis)
  249. )
  251. ### Run Batch Correction
  253. #Pre-allocation
  254. flattenedData = c()
  255. my_lengths = list()
  257. #For each subject, flatten out the vector to prepare for the batch correction
  258. for (sub in 1:length(filteredPheno$subid)){
  259. currVec = c()
  260. for (m in 1:dim_output){
  261. currVec = c(currVec, as.vector(unlist(turbData[[sub]][m])))
  262. }
  263. my_lengths[sub] = length(currVec)
  264. flattenedData = cbind(flattenedData, currVec)
  265. }
  267. #Vector-length of flattened participant data
  268. length_flattened_data = dim(flattenedData)[1]
  269. #Remove Lambdas from flattened data for the batch correction. Add the lambda data later.
  270. lambda_flattened = flattenedData[(length_flattened_data-nLambda+1):length_flattened_data,][,1]
  271. flattenedData = flattenedData[1:(length_flattened_data-nLambda),]
  273. #Batch variable
  274. batch = filteredPheno$site
  275. #Running the batch correction with the ComBat function
  276. model = model.matrix(~ mean_fd + sex + age_years + diagnosis + diagnosis:age_years, data = filteredPheno)
  277. correctedData = ComBat(dat=flattenedData, mod=model, batch=batch, par.prior=TRUE)
  279. #Testing if Batch correction worked with PCA
  280. #Extracting first PC from uncorrected data
  281. pca_uncorrected = prcomp(t(flattenedData), center = TRUE, scale. = TRUE)
  282. first_PC_uncorrected = pca_uncorrected$x[,1]
  284. #Run Anova for uncorrected PCA
  285. anova_model = aov(first_PC_uncorrected ~ batch)
  286. summary(anova_model)
  287. ```
  288. highly significant (p < 0.001)
  290. ```{r}
  291. #Extracting first PC from corrected data
  292. pca_corrected = prcomp(t(correctedData), center = TRUE, scale. = TRUE)
  293. first_PC_corrected = pca_corrected$x[,1]
  295. #Run Anova for corrected PCA
  296. anova_model = aov(first_PC_corrected ~ batch)
  297. summary(anova_model)
  298. ```
  299. Batch correction worked (p = 0.395)
  300. Plotting the results of the Batch correction
  302. ```{r}
  303. #Plotting results of the PCA before Batch correction
  304. PC_data = data.frame(batch, first_PC_corrected, first_PC_uncorrected)
  305. PC_beforeBatchCorrection = ggplot(PC_data, aes(x = batch, y = first_PC_uncorrected, color = batch)) +
  306. geom_scatterbox() +
  307. ylab("PC1") +
  308. xlab("Site") +
  309. guides(colour = "none") +
  310. easy_rotate_x_labels(angle = 45, side = "right") +
  311. ggtitle("Before batch correction") +
  312. theme(
  313. axis.title = element_text(size = 24), # axis labels
  314. axis.text = element_text(size = 18), # x and y axis tick labels
  315. plot.title = element_text(size = 24, hjust = 0.5) # title centered and larger
  316. )
  318. PC_beforeBatchCorrection
  319. ##ggsave(paste0(figures_directory,"/supplementary/PC_before_batch.pdf"), PC_beforeBatchCorrection, width = 16, height = 8)
  320. ```
  322. ```{r}
  323. #And after the batch correction
  324. PC_afterBatchCorrection = ggplot(PC_data, aes(x = batch, y = first_PC_corrected, color = batch)) +
  325. geom_scatterbox() +
  326. ylab("PC1") +
  327. xlab("Site") +
  328. guides(colour = "none") +
  329. easy_rotate_x_labels(angle = 45, side = "right") +
  330. ggtitle("After batch correction") +
  331. theme(
  332. axis.title = element_text(size = 24), # axis labels
  333. axis.text = element_text(size = 18), # x and y axis tick labels
  334. plot.title = element_text(size = 24, hjust = 0.5) # title centered and larger
  335. )
  337. PC_afterBatchCorrection
  338. ##ggsave(paste0(figures_directory,"/supplementary/PC_after_batch.pdf"), PC_afterBatchCorrection, width = 16, height = 8)
  339. ```
  341. ```{r}
  342. ### Reshaping
  343. # Reshaping the batch-corrected data into the different measures
  345. # Pre-allocation
  346. correctedTurbData = list()
  347. turbDataTest = list() #list to test if the dimensions line up with the initial data
  348. intermediate_num = (dim_node+nLambda+dim_RSN)
  350. # Loop over all subjects to save the batch-corrected data
  351. for (k in 1:length(filteredPheno$subid)){
  352. Turbulence = correctedData[1:nLambda,k]
  353. TurbulenceTest = flattenedData[1:nLambda,k]
  355. TurbulenceNode = matrix(correctedData[(nLambda+1):(dim_node+nLambda),k],nrow = nParcels, ncol = nLambda)
  356. TurbulenceNodeTest = matrix(flattenedData[(nLambda+1):(dim_node+nLambda),k],nrow = nParcels, ncol = nLambda)
  358. TurbulenceRSN = matrix(correctedData[(dim_node+nLambda+1):(dim_node+nLambda+dim_RSN),k], nrow = dim_RSN/nLambda, ncol = nLambda)
  359. TurbulenceRSNTest = matrix(flattenedData[(dim_node+nLambda+1):(dim_node+nLambda+dim_RSN),k], nrow = dim_RSN/nLambda, ncol = nLambda)
  361. InformationTransfer = correctedData[(intermediate_num+1):(intermediate_num+nLambda),k]
  362. InformationTransferTest = flattenedData[(intermediate_num+1):(intermediate_num+nLambda),k]
  364. InformationCascade = correctedData[intermediate_num+nLambda+1,k]
  365. InformationCascadeTest = flattenedData[intermediate_num+nLambda+1,k]
  367. InformationCascadeflow = correctedData[(intermediate_num+nLambda+2):(intermediate_num+(2*(nLambda))),k]
  368. InformationCascadeflowTest = flattenedData[(intermediate_num+nLambda+2):(intermediate_num+(2*(nLambda))),k]
  370. GlobalKuramoto = correctedData[intermediate_num+(2*(nLambda)+1),k]
  371. GlobalKuramotoTest = flattenedData[intermediate_num+(2*(nLambda)+1),k]
  373. Metastability = correctedData[intermediate_num+(2*(nLambda)+2),k]
  374. MetastabilityTest = flattenedData[intermediate_num+(2*(nLambda)+2),k]
  376. mean_KLOP = matrix(correctedData[(intermediate_num+(2*(nLambda)+3)):(dim(correctedData)[1]),k], ncol = nParcels, nrow = nLambda)
  377. mean_KLOPTest = matrix(flattenedData[(intermediate_num+(2*(nLambda)+3)):(dim(correctedData)[1]),k], ncol = nParcels, nrow = nLambda)
  379. Lambda = lambda_flattened
  380. LambdaTest = lambda_flattened
  382. subjectData = list(Turbulence,TurbulenceNode,TurbulenceRSN,InformationTransfer,InformationCascade,InformationCascadeflow,GlobalKuramoto,Metastability,mean_KLOP, Lambda)
  383. subjectDataTest = list(TurbulenceTest,TurbulenceNodeTest,TurbulenceRSNTest,InformationTransferTest,InformationCascadeTest,InformationCascadeflowTest,GlobalKuramotoTest,MetastabilityTest, mean_KLOPTest,LambdaTest)
  385. currentID = filteredPheno$subid[k]
  386. correctedTurbData[[as.character(currentID)]] = subjectData
  387. turbDataTest[[as.character(currentID)]] = subjectDataTest
  388. }
  390. #Extracting each measures
  391. amplTurb = t(as.data.frame(lapply(correctedTurbData, function(row) row[[1]])))
  392. nodeTurb = lapply(correctedTurbData, function(row) row[[2]])
  393. rsnTurb = lapply(correctedTurbData, function(row) row[[3]])
  394. transfer = t(as.data.frame(lapply(correctedTurbData, function(row) row[[4]])))
  395. cascade = as.numeric(unlist(lapply(correctedTurbData, function(row) row[[5]])))
  396. cascadeflow = t(as.data.frame(lapply(correctedTurbData, function(row) row[[6]])))
  397. gloKur = as.numeric(unlist(lapply(correctedTurbData, function(row) row[[7]])))
  398. meta = as.numeric(unlist(lapply(correctedTurbData, function(row) row[[8]])))
  399. mean_KLOP = t(as.data.frame(lapply(correctedTurbData, function(row) apply(row[[9]],1, FUN = mean))))
  400. nodeSync = lapply(correctedTurbData, function(row) t(row[[9]]))
  401. subid = as.factor(names(correctedTurbData))
  403. #Defining global datameasures
  404. global_data = data.frame(
  405. global_kuramoto = gloKur,
  406. metastability = meta,
  407. cascade = cascade,
  408. filteredPheno
  409. )
  411. ### Transforming ADOS values into percentage scores ###
  412. # Adding ADOS columns to the data frame
  413. global_data = global_data %>%
  414. mutate(
  415. ados2_socaff_css = NaN,
  416. ados2_rrb_css = NaN,
  417. ados2_total_css = NaN,
  418. ados_comb = NaN,
  419. ados_version = NaN
  420. )
  422. #Loop over subjects to compute the ADOS percentage score
  423. for (i in 1:nrow(global_data)){
  424. if (!is.na(global_data$ados2_total[i]) & !is.na(global_data$ados_module[i])){
  425. if (global_data$ados_module[i] == 4){
  426. denominator = 29
  427. }
  428. else{
  429. denominator = 28}
  430. global_data$ados_comb[i] = (global_data$ados2_total[i]/denominator)*100
  431. global_data$ados_version[i] = "ados"
  432. }
  433. else if(!is.na(global_data$ados_total[i]) & !is.na(global_data$ados_module[i])){
  434. denominator = 32
  435. global_data$ados_comb[i] = (global_data$ados_total[i]/denominator)*100
  436. global_data$ados_version[i] = "ados"}
  437. }
  439. # Filtering for subjects with existing total score
  440. raw_ados = global_data %>% dplyr::filter(!is.na(ados_comb))
  441. raw_ados$ados_version = as.factor(raw_ados$ados_version)
  442. ```
  444. Demographic data
  446. ```{r}
  447. ASD_data = filteredPheno %>% filter(diagnosis == "Autism")
  448. TD_data = filteredPheno %>% filter(diagnosis == "Control")
  450. #Ranges
  451. demo_summary = filteredPheno %>%
  452. filter(diagnosis %in% c("Autism", "Control")) %>%
  453. group_by(diagnosis) %>%
  454. summarise(
  455. age_min = min(age_years, na.rm = TRUE),
  456. age_max = max(age_years, na.rm = TRUE),
  457. age_mean = mean(age_years, na.rm = TRUE),
  458. age_sd = sd(age_years, na.rm = TRUE),
  460. FIQ_min = min(fiq, na.rm = TRUE),
  461. FIQ_max = max(fiq, na.rm = TRUE),
  462. FIQ_mean = mean(fiq, na.rm = TRUE),
  463. FIQ_sd = sd(fiq, na.rm = TRUE),
  465. VIQ_min = min(viq, na.rm = TRUE),
  466. VIQ_max = max(viq, na.rm = TRUE),
  467. VIQ_mean = mean(viq, na.rm = TRUE),
  468. VIQ_sd = sd(viq, na.rm = TRUE),
  470. PIQ_min = min(piq, na.rm = TRUE),
  471. PIQ_max = max(piq, na.rm = TRUE),
  472. PIQ_mean = mean(piq, na.rm = TRUE),
  473. PIQ_sd = sd(piq, na.rm = TRUE)
  474. )
  477. ### --------- age --------- ###
  478. #Histograms
  479. #ASD
  480. asd_age_hist = ggplot(filteredPheno[filteredPheno$diagnosis == "Autism",], aes(x = age_years)) +
  481. geom_histogram(alpha = 0.8, position = "identity", bins = 70, fill = "#ff7f0e", color = "black") +
  482. ylim(c(0,50))+
  483. xlim(c(0,65)) +
  484. theme_minimal() +
  485. theme(
  486. axis.title = element_blank(),
  487. axis.text = element_blank()
  488. )
  490. #TD
  491. td_age_hist = ggplot(filteredPheno[filteredPheno$diagnosis == "Control",], aes(x = age_years)) +
  492. geom_histogram(alpha = 0.8, position = "identity", bins = 70, fill = "#1f77b4", color = "black") +
  493. ylim(c(0,50))+
  494. xlim(c(0,65)) +
  495. theme_minimal() +
  496. theme(
  497. axis.title = element_blank(),
  498. axis.text = element_blank()
  499. )
  501. #Both
  502. comb_age_hist = ggplot(filteredPheno, aes(x = age_years, fill = diagnosis)) +
  503. geom_histogram(alpha = 0.5, position = "identity", bins = 70, color = "black") +
  504. scale_fill_manual(values = c("Autism" = "#ff7f0e", "Control" = "#1f77b4")) +
  505. ylim(c(0,50))+
  506. xlim(c(0,65)) +
  507. theme_minimal() +
  508. theme(
  509. axis.title = element_blank(),
  510. axis.text = element_blank(),
  511. legend.position = "none"
  512. )
  514. ### --- Printing and Saving -- ###
  515. print(td_age_hist)
  516. print(asd_age_hist)
  517. print(comb_age_hist)
  518. #ggsave(paste0(figures_directory, "/supplementary/age/td_age_hist.jpg"), width = 10, height = 6, td_age_hist, dpi = 300)
  519. #ggsave(paste0(figures_directory, "/supplementary/age/asd_age_hist.jpg"),width = 10, height = 6 ,asd_age_hist, dpi = 300)
  520. #ggsave(paste0(figures_directory, "/supplementary/age/comb_age_hist.jpg"),width = 10, height = 6, comb_age_hist, dpi = 300)
  522. ```
  526. Overview of the Turbulence measures in the data list (with dimensions in brackets):
  527. 1: Amplitude Turbulence [10,1], 2: Parcel-level Turbulence [1054,10], 3: RSN Turbulence [8,10]
  528. 4: Information Transfer [10,1], 5: Information Cascade [1], 6: Information Cascadeflow [9,1]
  529. 7: GolbalKuramoto [1], 8: Metastability [1], 9: Node_Synchronization, 10: Lambda values [1,10]
  531. Lambda values range from high to low, so small to large distances (1=0.28, 2=0.25, ... ,10 = 0.01)
  533. 1.0: Statistical analysis of whole-brain summary statistics: Global Synchronization, Metastability
  534. ```{r, warning=FALSE, message=FALSE}
  535. #Variable handling and pre-allocation
  536. counter = 0
  537. global_measures = c("global_kuramoto", "metastability")
  538. axes_labels = c("Global Synchronization", "Metastability")
  539. global_statistics = list()
  541. #Loop over all global measures
  542. for (curr_measure in global_measures){
  543. counter = counter+1
  545. ### TESTS
  547. # 1) Diagnosis as IV
  548. # ANOVAs
  549. diag_mod = lm(formula = global_data[[curr_measure]] ~ mean_fd + sex + age_years + diagnosis + diagnosis:age_years, data = global_data) #create model
  550. diag_ANOVA = anova(diag_mod)
  551. #Compute Cohens D:
  552. #Group sizes
  553. n1 = sum(global_data$diagnosis == "Control")
  554. n2 = sum(global_data$diagnosis == "Autism")
  555. diag_f = diag_ANOVA["diagnosis", "F value"]
  556. diag_d = round(sqrt(diag_f * (n1 + n2) / (n1 * n2)), 3)
  557. diag_eta = c(round(effectsize(diag_ANOVA)$Eta2_partial,3),NA)
  559. # 2) ADOS as IV
  560. #ANOVAs
  561. ados_mod = lm(formula = raw_ados[[curr_measure]] ~ mean_fd + sex + age_years + ados_comb + ados_comb:age_years,
  562. data = raw_ados)
  563. ados_ANOVA = anova(ados_mod)
  564. ados_f = ados_ANOVA["ados_comb", "F value"]
  565. ados_d = round(sqrt(ados_f * (n1 + n2) / (n1 * n2)), 3)
  567. # Correlations of ADOS with measures
  568. # Compute correlation test
  569. curr_corr = cor.test(raw_ados[[curr_measure]], raw_ados$ados_comb)
  571. # Correlations of FD with measures differences between ASD and TD
  572. curr_ASD = global_data %>% filter(diagnosis == "Autism")
  573. curr_TD = global_data %>% filter(diagnosis == "Control")
  574. ASD_corr = cor.test(curr_ASD[[curr_measure]], curr_ASD$mean_fd)
  575. TD_corr = cor.test(curr_TD[[curr_measure]], curr_TD$mean_fd)
  576. curr_paired_r = paired.r(ASD_corr$estimate, TD_corr$estimate, NULL, nrow(curr_ASD), nrow(curr_TD))
  577. curr_paired_p = curr_paired_r$p
  578. curr_paired_z = curr_paired_r$z
  579. curr_paired_stats = data.frame(
  580. ASD_corr$estimate, TD_corr$estimate, curr_paired_z, curr_paired_p)
  582. #Merging output
  583. global_statistics[[curr_measure]] = list(data.frame(diag_ANOVA),
  584. paste("d value diagnosis",diag_d),
  585. data.frame(ados_ANOVA),
  586. paste("d value ados",ados_d),
  587. data.frame(curr_corr$estimate,
  588. curr_corr$p.value,
  589. curr_corr$statistic),
  590. curr_paired_stats)
  592. ### PLOTS ###
  593. # Boxplot of measure dx diagnosis
  594. curr_boxplot = ggplot(data = global_data, aes(x = diagnosis, y = global_data[[curr_measure]], colour = diagnosis)) +
  595. geom_scatterbox() +
  596. ylab(axes_labels[counter]) +
  597. xlab("Diagnosis") +
  598. guides(colour = "none") +
  599. scale_x_discrete(labels = c("Autism" = "Autism", "Control" = "TD"))+
  600. theme_minimal() +
  601. ylim(NA, max(global_data[[curr_measure]])+0.1) +
  602. scale_color_manual(values = c("Control" = "#1f77b4", "Autism" = "#ff7f0e")) +
  603. theme(
  604. axis.title = element_blank(),
  605. axis.text = element_blank(),
  606. )
  608. #Correlation of the measure dx ados
  609. curr_corr_plot = ggplot(raw_ados, aes(x = ados_comb, y = .data[[curr_measure]])) +
  610. geom_jitter(size = 5) +
  611. geom_smooth(method = "lm", linewidth = 3) +
  612. theme_minimal() +
  613. theme(
  614. axis.title = element_blank(),
  615. axis.text = element_blank(),
  616. )
  617. curr_corr_plot = curr_corr_plot +
  618. theme(plot.margin = margin(0, 15, 0, 0))
  620. # Printing and saving
  621. print(curr_boxplot)
  622. print(curr_corr_plot)
  623. ##ggsave(paste0(figures_directory,"/Figure_3/",curr_measure,".jpg"), curr_boxplot, width = 8, height = 6)
  624. #ggsave(paste0(figures_directory,"/supplementary/S3/",curr_measure,"_ADOS_corr.jpg"), curr_corr_plot, width = 10, height = 6, dpi = 300)
  625. ##ggsave(paste0(figures_directory,"/Figure_7/",curr_measure,".jpg"), curr_boxplot, width = 8, height = 6)
  626. }
  629. #Change the list names
  630. for (c in names(global_statistics)){
  631. names(global_statistics[[c]]) = c("ANOVA_diagnosis", "diagnosis_d", "ANOVA_ados", "ados_d", "correlation_ados")}
  632. ```
  634. 2.0 Statistical analysis of local measures: Scale-specific Synchronization, Amplitude Turbulence, Information Cascadeflow, Information Transfer
  636. ```{r, warning=FALSE, message=FALSE}
  637. #Variable handling and pre-allocation
  638. local_statistics = list()
  639. turbulence_ados_stats = list()
  640. corr_meanFD_summary = list()
  641. local_measures = c("turbulence", "transfer", "cascadeflow", "scale_synchronization")
  642. local_vectors = list(
  643. "turbulence" = unlist(lapply(amplTurb, function(row) row[1])),
  644. "transfer" = unlist(lapply(transfer, function(row) row[1])),
  645. "cascadeflow" = unlist(lapply(cascadeflow, function(row) row[1])),
  646. "scale_synchronization" = unlist(lapply(mean_KLOP, function(row) row[1]))
  647. )
  648. plot_labels = c("Amplitude turbulence", "Information transfer", "Information cascade flow", "Scale-specific synchronization")
  649. counter = 0
  650. curr_measure = "turbulence"
  651. for (curr_measure in local_measures){
  652. counter = counter +1
  653. if (curr_measure == "cascadeflow"){
  654. nLambda = 9
  655. Lambda = lambda_flattened[-1]
  656. } else{
  657. nLambda = 10
  658. Lambda = lambda_flattened
  659. }
  660. #Data initialization
  661. curr_data = data.frame(
  662. subid = as.factor(rep(subid, times = nLambda)),
  663. lambda = as.numeric(rep(Lambda, each = length(global_data$subid))),
  664. measure = local_vectors[[curr_measure]],
  665. diagnosis = factor(rep(global_data$diagnosis, times = nLambda)),
  666. sex = factor(rep(global_data$sex, times = nLambda)),
  667. mean_fd = rep(global_data$mean_fd, times = nLambda),
  668. age_years = rep(global_data$age_years, times = nLambda),
  669. ados = rep(global_data$ados_comb, times = nLambda))
  671. #1) Diagnosis as IV
  672. #ANOVA
  673. diag_mod = lmerTest::lmer(formula = measure ~ mean_fd + sex + age_years + lambda + diagnosis + diagnosis:lambda + diagnosis:age_years + (1+lambda|subid), data = curr_data) #Model initialization
  674. diag_ANOVA = anova(diag_mod, type = 1) #Step-wise Anova
  675. diag_eta = round(effectsize(diag_ANOVA)$Eta2_partial,3)
  676. diag_emm_curr = emmeans(diag_mod, ~ diagnosis | lambda, at = list(lambda = sort(unique(curr_data$lambda)))) #Posthoc tests for each Lambda level
  677. diag_stat = summary(pairs(diag_emm_curr))
  678. curr_lambda = diag_stat$lambda
  679. diag_p = diag_stat$p.value
  680. diag_d = summary(eff_size(diag_emm_curr, sigma = sigma(diag_mod), edf = df.residual(diag_mod))) #Cohens d in linear models is d = mean_1 - mean2 / sigma, where sigma is residual standard deviation
  681. diag_d = diag_d$effect.size
  682. diag_means = summary(diag_emm_curr)$emmean
  683. means_ASD = diag_means[seq(1,length(diag_p)*2, by = 2)]
  684. means_TD = diag_means[seq(2,length(diag_p)*2, by = 2)]
  685. diag_stat_summary = data.frame(
  686. means_ASD = means_ASD,
  687. means_TD = means_TD,
  688. means_diff = means_ASD - means_TD,
  689. p = p.adjust(diag_p, method = "BH"),
  690. d = diag_d,
  691. lambda = curr_lambda
  692. )
  694. #ANOVA results saving
  695. diag_ANOVA_results = cbind(data.frame(diag_ANOVA),
  696. eta = diag_eta)
  698. #2) ADOS as IV
  699. #Filtering data to only autisms with existing ados percentage score
  700. curr_ados_data = curr_data %>% dplyr::filter(!is.na(ados))
  701. ados_mod = lmerTest::lmer(formula = measure ~ mean_fd + sex + age_years + lambda + ados + ados:lambda + ados:age_years + (1+lambda|subid), data = curr_ados_data) #Model initialization
  702. ados_ANOVA = anova(ados_mod, type = 1) #Step-wise Anova
  703. ados_trends = emtrends(ados_mod, ~ lambda, var = "ados",
  704. at = list(lambda = sort(unique(curr_data$lambda))))
  705. ados_stats = summary(ados_trends, infer = TRUE)
  707. # Correlations of measures with ados
  708. curr_ados_data = curr_data %>%
  709. group_by(subid, age_years, diagnosis, ados) %>%
  710. summarise(overall_mean = mean(measure))
  711. curr_ados_data = curr_ados_data %>% filter(!is.na(ados))
  712. curr_ados_corr = cor.test(curr_ados_data$ados, curr_ados_data$overall_mean)
  714. # Saving in data frame
  715. local_statistics[[curr_measure]] = list(
  716. diag_stat_summary,
  717. diag_ANOVA_results,
  718. data.frame(ados_ANOVA),
  719. ados_stats,
  720. data.frame(ados_corr_p = curr_ados_corr$p.value,
  721. ados_corr_r = curr_ados_corr$estimate))
  723. # Investigating age effects
  724. curr_age_data = curr_data %>% group_by(subid, age_years, diagnosis) %>% summarise(overall_mean = mean(measure))
  725. curr_age_data_ASD = curr_age_data %>% filter(diagnosis == "Autism")
  726. curr_age_data_TD = curr_age_data %>% filter(diagnosis == "Control")
  727. ASD_cor_test = cor.test(curr_age_data_ASD$age_years, curr_age_data_ASD$overall_mean)
  728. TD_cor_test = cor.test(curr_age_data_TD$age_years, curr_age_data_TD$overall_mean)
  729. a = paired.r(ASD_cor_test$estimate,TD_cor_test$estimate, NULL, 504,505)
  731. #Investigating effects of meanFD on metrics
  732. curr_meanFD_data = curr_data %>%
  733. group_by(subid, mean_fd, diagnosis) %>%
  734. summarise(overall_mean = mean(measure))
  736. #Correlation meanFD ~ autism
  737. asd_data = curr_meanFD_data %>%
  738. filter(diagnosis == "Autism")
  739. corr_autism_meanFD = cor.test(asd_data$overall_mean, asd_data$mean_fd)
  741. #Correlation meanFD ~ TD
  742. control_data = curr_meanFD_data %>%
  743. filter(diagnosis == "Control")
  744. corr_td_meanFD = cor.test(control_data$overall_mean, control_data$mean_fd)
  746. #Fisher z-test
  747. curr_paired_r = paired.r(corr_autism_meanFD$estimate, corr_td_meanFD$estimate, NULL,
  748. nrow(asd_data), nrow(control_data))
  749. curr_paired_p = curr_paired_r$p
  750. curr_paired_z = curr_paired_r$z
  751. curr_paired_stats = data.frame(
  752. corr_autism_meanFD$estimate, corr_td_meanFD$estimate, curr_paired_z, curr_paired_p)
  754. #Saving
  755. corr_meanFD_summary[[curr_measure]] = curr_paired_stats
  757. ### PLOTS
  758. #Inverting Lambda-levels
  759. curr_data$lambda = as.factor(curr_data$lambda) #Only for the plots!
  760. #curr_data$lambda = factor(curr_data$lambda, levels = rev(levels(curr_data$lambda)))
  762. #Measure dx meanFD
  763. curr_meanFD_plot = ggplot(data = curr_meanFD_data, aes(x = mean_fd, y = overall_mean, colour = diagnosis, group = diagnosis)) +
  764. geom_jitter() +
  765. geom_smooth(method = "lm") +
  766. ylab(plot_labels[counter]) +
  767. xlab("mean framewise displacement") +
  768. scale_color_manual(values = c("Control" = "#1f77b4", "Autism" = "#ff7f0e")) +
  769. theme_minimal() +
  770. theme(
  771. axis.title = element_text(size = 30),
  772. axis.text = element_text(size = 30),
  773. legend.text = element_blank(),
  774. legend.title = element_blank(),
  775. legend.position = "none")
  777. #Boxplot of measure and Lambda by diagnosis
  778. curr_dx_Lambda_by_diagnosis_boxplot = ggplot(data = curr_data,
  779. aes(x = lambda,
  780. y = measure,
  781. color = diagnosis)) +
  782. geom_jitter(position = position_jitterdodge(jitter.width = 0.5, dodge.width = 1), alpha = 0.5, size = 1.5) +
  783. geom_boxplot(aes(group = interaction(lambda, diagnosis)), colour = "black",fill = NA,position=position_dodge(width=1),outlier.shape = NA, size = 1.5) +
  784. ylab(plot_labels[counter]) +
  785. xlab("Lambda (in 1/mm)") +
  786. scale_color_manual(values = c("Control" = "#1f77b4", "Autism" = "#ff7f0e")) +
  787. theme_minimal() +
  788. theme(
  789. axis.title = element_text(size = 30),
  790. axis.text = element_text(size = 30),
  791. legend.text = element_blank(),
  792. legend.title = element_blank(),
  793. legend.position = "none")
  795. #Measure dx age
  796. curr_age_plot = ggplot(data = curr_age_data, aes(x = age_years, y = overall_mean, colour = diagnosis, group = diagnosis)) +
  797. geom_jitter(size = 4) +
  798. geom_smooth(method = "lm", linewidth = 4) +
  799. scale_color_manual(values = c("Control" = "#1f77b4", "Autism" = "#ff7f0e")) +
  800. theme_minimal() +
  801. theme(
  802. axis.title = element_blank(),
  803. axis.text = element_blank(),
  804. legend.position = "none")
  806. #Measure dx ADOS
  807. curr_ados_plot = ggplot(data = curr_ados_data, aes(x = ados, y = overall_mean)) +
  808. geom_jitter(size = 5) +
  809. geom_smooth(method = "lm", linewidth = 3) +
  810. theme_minimal() +
  811. theme(
  812. axis.text = element_blank(),
  813. axis.title = element_blank(),
  814. legend.position = "none")
  816. # Correlation interaction for Turbulence
  817. if (curr_measure == "turbulence"){
  818. lambdas_to_plot = c(0.01,0.28)
  819. for (curr_lambda in lambdas_to_plot){
  820. curr_lam_data = curr_data %>% filter(lambda == curr_lambda)
  821. curr_ados_lambda_plot = ggplot(data = curr_lam_data, aes(x = ados, y = measure)) +
  822. geom_jitter(size = 5) +
  823. geom_smooth(method = "lm", linewidth = 3) +
  824. theme_minimal() +
  825. theme(
  826. axis.text = element_blank(),
  827. axis.title = element_blank(),
  828. legend.position = "none")
  830. #Correlation Tests
  831. # Correlations of measures with ados
  832. curr_ados_corr = cor.test(curr_lam_data$ados, curr_lam_data$measure)
  833. turbulence_ados_stats[[paste0(curr_lambda)]] = data.frame(
  834. ados_corr_p = curr_ados_corr$p.value,
  835. ados_corr_r = curr_ados_corr$estimate)
  836. print(curr_ados_lambda_plot)
  837. ##ggsave(paste0(figures_directory,"/supplementary/S3/",curr_measure,"_",curr_lambda,"_ados_corr.jpg"), curr_ados_lambda_plot, height = 6, width = 10, dpi = 300)
  838. }
  839. }
  841. #Printing and saving
  842. print(curr_dx_Lambda_by_diagnosis_boxplot)
  843. print(curr_age_plot)
  844. #print(curr_ados_plot)
  845. #print(curr_meanFD_plot)
  846. ##ggsave(paste0(figures_directory,"/Figure_6/box_",curr_measure,".pdf"), curr_dx_Lambda_by_diagnosis_boxplot, height = 12, width = 16)
  847. ##ggsave(paste0(figures_directory,"/Figure_6/diff_",curr_measure,".pdf"), curr_diff_plot, height = 12, width = 16)
  848. ##ggsave(paste0(figures_directory,"/supplementary/S3/",curr_measure,"_ados_corr.jpg"), curr_ados_plot, height = 6, width = 10, dpi = 300)
  849. #ggsave(paste0(figures_directory,"/supplementary/S5/",curr_measure,"_age_corr.jpg"), curr_age_plot, height = 6, width = 10, dpi = 300)
  850. }
  851. ```
  854. 3. Network-wise analysis
  856. In this part of the analysis, we want to investigate if there are Turbulence patterns of ASD vs. TD in different RSNs.
  857. To do so, we leverage the node-turbulence and RSN turbulence matrices.
  859. 3.1. Node-wise RSN analysis. Starting from individual node turbulence matrices (nNodes*nLambda). Computing average for ASD and TD and calculate the difference.
  860. The highest 10% values of differences are extracted for each lambda. These 100 nodes are then associated to the RSNs.
  862. ```{R, message=FALSE, warning=FALSE}
  863. ### Loading and initializing data
  864. #Mapping from RSN number to label:
  865. rsns = c("VIS", "SOM", "DAN", "VAN", "LIM", "CON", "DMN", "SUB")
  866. rsnStats = list()
  867. rsn_labels = unlist(readMat(file.path(root, "code/labels.mat"))) #Loading RSN-labels of parcellation
  868. schaef_regions = read.table(file.path(root, "code/Schaefer2018_1000Parcels_7Networks_order.txt")) #Load Schaefer region names
  869. tian_regions = read.table(file.path(root, "code/Tian_Subcortex_S4_3T_label.txt")) #Load Tian region names
  870. Lambda = c(0.28,0.25,0.22,0.19,0.16,0.13,0.10,0.07,0.04,0.01)
  872. #Pre-allocation
  873. parcel_indices = list()
  874. parcel_ind_TD_higher = list()
  875. nodeTurb_residuals = array(NA,dim = c(nParcels, nLambda, nrow(filteredPheno)))
  877. # Testing if parcel level turbulence is associated to FD differently between ASD and TD
  878. Autism_data = nodeTurb[which(filteredPheno$diagnosis == "Autism")]
  879. TD_data = nodeTurb[which(filteredPheno$diagnosis == "Control")]
  880. my_ps = data.frame()
  881. my_zs = data.frame()
  882. for (curr_parcel in 1:nParcels){
  883. for (curr_lambda in 1:nLambda){
  885. ### -------- Correlations of parcel-level turbulence with FD -------- ###
  887. #Correlation of parcel turbulence with FD for ASD
  888. currASD_data = drop_na(data.frame(
  889. brain = sapply(Autism_data, function(x) {x[curr_parcel,curr_lambda]}),
  890. fd = filteredPheno[filteredPheno$diagnosis == "Autism",]$mean_fd))
  891. curr_cor_ASD = cor.test(currASD_data$brain,currASD_data$fd)
  893. #Correlation of parcel turbulence with FD for TD
  894. currTD_data = drop_na(data.frame(
  895. brain = sapply(TD_data, function(x) {x[curr_parcel,curr_lambda]}),
  896. fd = filteredPheno[filteredPheno$diagnosis == "Control",]$mean_fd))
  897. curr_cor_TD = cor.test(currTD_data$brain,currTD_data$fd)
  899. # Ztest of difference
  900. curr_paired = paired.r(curr_cor_ASD$estimate, curr_cor_TD$estimate, NULL, nrow(currASD_data), nrow(currTD_data))
  901. my_ps[curr_parcel,curr_lambda] = curr_paired$p
  902. my_zs[curr_parcel,curr_lambda] = curr_paired$z
  904. ### ------ computation of parcel-level turbulence residuals after covariate regression ------ ###
  905. curr_Parcel_Turb_data = drop_na(data.frame(
  906. parcel_turb = sapply(nodeTurb, function(x) {x[curr_parcel,curr_lambda]}),
  907. fd = filteredPheno$mean_fd,
  908. sex = as.factor(filteredPheno$sex),
  909. age = filteredPheno$age_years))
  911. my_model = lm(formula = parcel_turb ~ fd + sex + age, data = curr_Parcel_Turb_data)
  912. curr_residuals = resid(my_model)
  913. nodeTurb_residuals[curr_parcel, curr_lambda,] = curr_residuals
  914. }
  915. }
  918. #Compute FDR
  919. my_ps_fdr = data.frame(matrix(
  920. p.adjust(flatten(my_ps), method = "BH"),
  921. nrow = nParcels,
  922. ncol = length(lambda_flattened)
  923. ))
  925. #write_csv(my_zs, paste0(root,"/code/z_values_meanFD_parcelTurb.csv"))
  926. #write_csv(my_ps_fdr, paste0(root,"/code/p_values_meanFD_parcelTurb.csv"))
  928. ### No significant differences in correlations between groups
  930. #Set data to node Turbulence residuals
  931. nSubjects = dim(nodeTurb_residuals)[3]
  932. nodeTurb_residuals = lapply(
  933. 1:nSubjects,
  934. function(s) nodeTurb_residuals[ , , s]
  935. )
  936. names(nodeTurb_residuals) = names(nodeTurb)
  938. #Extracting masks of ASD and control
  939. ASD_labels = which(filteredPheno$diagnosis == "Autism")
  940. TD_labels = which(filteredPheno$diagnosis == "Control")
  941. ASD_node = nodeTurb_residuals[ASD_labels]
  942. TD_node = nodeTurb_residuals[TD_labels]
  943. parcel_labels = c(unlist(schaef_regions$V2), tian_regions$V1)
  945. #Computing means and effect sizes
  946. ASD_mean = Reduce("+", ASD_node) / length(ASD_labels)
  947. TD_mean = Reduce("+", TD_node) / length(TD_labels)
  948. ASD_3d = simplify2array(ASD_node)
  949. TD_3d = simplify2array(TD_node)
  950. n_asd = length(ASD_node)
  951. n_td = length(TD_node)
  953. #Standard deviation
  954. pooled_sd = sqrt((
  955. apply(ASD_3d, c(1,2), var)*(n_asd-1) +
  956. apply(TD_3d, c(1,2), var)*(n_td-1)
  957. )/(n_asd+n_td-2))
  959. #Save Cohen's D
  960. turb_CohensD = (TD_mean-ASD_mean)/pooled_sd
  962. ### Network-wise analysis and radar plots
  963. #Absolute node-diff
  964. node_diff = abs(TD_mean-ASD_mean)
  966. # Number of parcels to be investigated
  967. n_Top_Parcels = 200
  969. #Extracting top 10% of absolute node-difference for each lamba
  970. ind = apply(node_diff, 2, function(col) {
  971. top = sort(col, decreasing = TRUE)[1:n_Top_Parcels]
  972. col %in% top})
  974. #Save in list
  975. parcel_indices = ind
  977. rsn_ind = apply(ind, 2, function(col) rsn_labels[col]) #Assign to each of these parcels its RSN according to the Schaefer parcellation
  978. counts = apply(rsn_ind, 2, table) #Compute for each Lambda its counts
  980. #Add missing values in the counts with 0:
  981. list_names = sapply(counts, names)
  982. for (i in 1:length(counts)){
  983. for (ii in 1:length(rsns)){
  984. if (!(as.character(ii) %in% list_names[[i]])){
  985. counts[[i]][as.character(ii)] = 0
  986. }
  987. }
  988. to_order = counts[[i]]
  989. counts[[i]] = to_order[order(as.numeric(names(to_order)))]
  990. }
  992. #Inverting the list to save in a data-frame.
  993. #Also, scale each count with the maximum number of parcels of that RSN. Result is a percentage score for each RSN
  994. n_parcels_RSN = as.numeric(table(rsn_labels))
  995. counts_new = matrix(nrow = length(Lambda), ncol = length(rsns))
  996. for (i in 1:length(counts)){
  997. counts_new[i,] = as.numeric(counts[[i]])#/n_parcels_RSN
  998. }
  999. #counts_new = counts_new*(1/max(counts_new))
  1000. counts = data.frame(Lambda, counts_new)
  1001. colnames(counts) = c("Lambda",rsns)
  1002. rownames(counts) = Lambda
  1004. ##### Directed analysis, i.e., checking which networks are higher in ASD vs. TD
  1006. #Compute directed node-difference
  1007. node_direct = TD_mean-ASD_mean
  1009. #Checking for each lambda, how many of the values are higher
  1010. #preallocation
  1011. ind_TD_higher = ind
  1012. rsn_diff = 0
  1013. for (lam in 1:10){
  1014. count = 0
  1015. curr_vec = node_direct[,lam]
  1016. for (parcel in 1:length(curr_vec)){
  1017. ind_TD_higher[parcel,lam] = FALSE
  1018. if (curr_vec[parcel] > 0 && ind[parcel,lam]){
  1019. count = count+1
  1020. ind_TD_higher[parcel,lam] = TRUE
  1021. }
  1022. }
  1023. rsn_diff[lam] = count
  1024. }
  1026. parcel_ind_TD_higher = ind_TD_higher
  1028. ### Plotting
  1030. color_vec = colorRampPalette(c("#1f77b4", "black", "#ff7f0e"))(10)
  1032. #1) Combined radarplot
  1033. combined_radar = ggradar(counts,
  1034. axis.label.size = 5,
  1035. legend.position = "right",
  1036. group.point.size = 4,
  1037. group.colours = (color_vec),
  1038. group.line.width = 1.2,
  1039. values.radar = c(0,max(counts)/2,max(counts)),
  1040. base.size = 10,
  1041. label.gridline.min = FALSE,
  1042. grid.mid = max(counts)/2,
  1043. grid.max = max(counts),
  1044. legend.text.size = 10,
  1045. legend.title = "Lambda",
  1046. background.circle.colour = "white",
  1047. background.circle.transparency = 0.2) +
  1048. theme(
  1049. #axis.title = element_text(size = 18), # axis labels
  1050. #axis.text = element_text(size = 18), # x and y axis tick labels
  1051. #plot.title = element_text(size = 18, hjust = 0.5), # title centered and larger
  1052. legend.text = element_blank(), # legend label font size
  1053. legend.title = element_blank(),
  1054. legend.position = "none")
  1056. #2) Barplot which shows the distribution of contribution of autism and TD as a function of lambda
  1057. # Adjust rsn_diff so the reference point is 50
  1058. dist_data = data.frame(
  1059. Lambda = rep(Lambda,2),
  1060. values = c((n_Top_Parcels - rsn_diff),rsn_diff),
  1061. cohort = rep(c("TD>Autism", "Autism>TD"), each = nLambda)
  1062. )
  1064. distribution_plot = ggplot(dist_data, aes(x = Lambda, y = values, fill = cohort)) +
  1065. geom_bar(stat = "identity") +
  1066. scale_fill_manual(values = c("TD>Autism" = "#00CED1", "Autism>TD" = "#A0522D")) +
  1067. scale_x_reverse() +
  1068. labs(
  1069. x = "Lambda (in 1/mm)",
  1070. y = "#parcels TD > ASD (in %)",
  1071. fill = "Distribution"
  1072. ) +
  1073. theme_minimal() +
  1074. theme(
  1075. axis.title = element_blank(),
  1076. axis.text = element_blank(),
  1077. axis.ticks.x = element_blank(),
  1078. axis.ticks.y = element_blank(),
  1079. plot.title = element_blank(),
  1080. legend.position = "none",
  1081. panel.background = element_rect(fill = "white", color = NA),
  1082. panel.grid.major.y = element_line(color = "grey75", size = 0.5),
  1083. panel.grid.minor.y = element_line(color = "grey75", size = 0.25),
  1084. panel.grid.minor.x = element_blank(),
  1085. panel.grid.major.x = element_blank()
  1086. )
  1089. #Printing and saving
  1090. print(distribution_plot)
  1091. print(combined_radar)
  1093. #ggsave(paste0(figures_directory,"/Figure_4/RSN_distribution_",n_Top_Parcels, "parcels.jpg"), distribution_plot, height = 12, width = 16, dpi = 300)
  1094. #ggsave(paste0(figures_directory,"/Figure_4/RSN_radar_",n_Top_Parcels, "parcels.jpg"), combined_radar, height = 6, width = 6,dpi = 300)
  1095. ```
  1097. 3.2 Association of 100 Nodes with biggest differences to Sydnor gradients
  1098. ```{R}
  1099. #First, we need to compute the rankings of the 1054 parcels. We do this with a python script. The input to the python script are the parcel-level effect sizes of both, the synchronization strength, as well as the turbulence. Here we save these into the directory in which the code of the python script operates.
  1101. #Turbulence data
  1102. currAutism_data = nodeTurb_residuals[which(filteredPheno$diagnosis == "Autism")]
  1103. currAutism_data = data.frame(Reduce("+", currAutism_data) / length(currAutism_data))
  1104. currTD_data = nodeTurb_residuals[which(filteredPheno$diagnosis == "Control")]
  1105. currTD_data = data.frame(Reduce("+", currTD_data) / length(currTD_data))
  1106. curr_d = turb_CohensD
  1108. #Save csv files
  1109. write_csv(data.frame(currAutism_data), paste0(root,"/data/gradient_analysis/", paste0("parcel_turb","_autism_1054", ".csv")))
  1110. write_csv(data.frame(currTD_data), paste0(root,"/data/gradient_analysis/", paste0("parcel_turb","_TD_1054", ".csv")))
  1111. write_csv(data.frame(curr_d), paste0(root,"/data/gradient_analysis/", paste0("parcel_turb","_CohensD_1054", ".csv")))
  1113. #Run the python script with the reticulate option
  1114. use_virtualenv("~/pyenvs/reticulate_arm", required = TRUE)
  1115. use_python("/Users/bened/pyenvs/reticulate_arm/bin/python", required = TRUE)
  1116. py_run_file(paste0(root, "/code/", "gradient_computations.py"))
  1117. ```
  1118. Association of parcel-level turbulence levels with S-A and functional gradients
  1120. ```{R}
  1121. #Loading in data
  1122. weighted_sydnor = read.csv(paste0(root, "/data/gradient_analysis/gradients_1054.csv"))
  1123. weighted_sydnor = weighted_sydnor[1:1000,c(1:11)]#only cortical regions
  1124. names(weighted_sydnor) = c("Anatomical T1w/T2w", "Functional Gradient", "Cortical Expansion", "Allometric Scaling",
  1125. "Aerobic Glycolysis", "CBF", "Gene Expression", "NeuroSynth", "Externopyramidization", "Cortical Thickness",
  1126. "Sensorimotor-Association Gradient")
  1127. plotnames = c("Anatomical_T1wT2w", "Functional_Gradient", "Cortical_Expansion", "Allometric_Scaling",
  1128. "Aerobic_Glycolysis", "CBF", "Gene_Expression", "NeuroSynth", "Externopyramidization", "Cortical_Thickness",
  1129. "Sensorimotor_Association_Gradient")
  1131. #Directions of the gradient for each measure: True from high to low, False from Low to high
  1132. gradient_directions = data.frame(
  1133. gradient = names(weighted_sydnor),
  1134. direction = c(FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE, TRUE)
  1135. )
  1137. #Since we are doing the comparison of the parcel-level turbulence with the S-A axes in the Schaefer1000 parcellation, we have to compute the parcels with the biggest differences again on that parcellation
  1138. rsn_labels = rsn_labels[1:1000] #Change labels to only the first 1000
  1139. parcel_labels = c(unlist(schaef_regions$V2))
  1140. node_cortical = lapply(nodeTurb_residuals, function(row) row[1:1000,])
  1141. ASD_node = node_cortical[ASD_labels]
  1142. TD_node = node_cortical[TD_labels]
  1144. #Computing means and effect sizes
  1145. ASD_mean = Reduce("+", ASD_node) / length(ASD_labels)
  1146. TD_mean = Reduce("+", TD_node) / length(TD_labels)
  1147. ASD_3d = simplify2array(ASD_node)
  1148. TD_3d = simplify2array(TD_node)
  1149. n_asd = length(ASD_node)
  1150. n_td = length(TD_node)
  1152. pooled_sd = sqrt((
  1153. apply(ASD_3d, c(1,2), var)*(n_asd-1) +
  1154. apply(TD_3d, c(1,2), var)*(n_td-1)
  1155. )/(n_asd+n_td-2))
  1157. #Save Cohen's D
  1158. turb_CohensD = (TD_mean-ASD_mean)/pooled_sd
  1160. ### Network-wise analysis and radar plots
  1161. #Absolute node-diff
  1162. node_diff = abs(TD_mean-ASD_mean)
  1164. #Extracting top 10% of absolute node-difference for each lamba
  1165. ind = apply(node_diff, 2, function(col) {
  1166. top = sort(col, decreasing = TRUE)[1:100]
  1167. col %in% top})
  1169. parcel_indices = ind
  1171. ### Directed analysis
  1172. #Compute directed node-difference
  1173. node_direct = TD_mean-ASD_mean
  1175. #Checking for each lambda, how many of the values are higher
  1176. #preallocation
  1177. ind_TD_higher = ind
  1178. rsn_diff = 0
  1179. for (lam in 1:10){
  1180. count = 0
  1181. curr_vec = node_direct[,lam]
  1182. for (parcel in 1:length(curr_vec)){
  1183. ind_TD_higher[parcel,lam] = FALSE
  1184. if (curr_vec[parcel] > 0 && ind[parcel,lam]){
  1185. count = count+1
  1186. ind_TD_higher[parcel,lam] = TRUE
  1187. }
  1188. }
  1189. rsn_diff[lam] = count
  1190. }
  1192. parcel_ind_TD_higher = ind_TD_higher
  1194. #### Assign to each of the 100 highest nodes a label, if TD is higher -> True
  1195. my_parcels = matrix(sapply(1:ncol(ind), function(j) {
  1196. parcel_ind_TD_higher[parcel_indices[, j], j]
  1197. }), nrow = 100, ncol = 10)
  1199. gradient_mod_comp = list()
  1201. ### Statistical analysis. Looping through the (selected) gradients to statistically test which model fits best
  1202. for (metric in c(2,11)){
  1203. #Pre-allocation
  1204. curr_csv = weighted_sydnor
  1205. data_sensory_assoc = data.frame()
  1206. #Loading data for that specific gradient
  1207. for (i in 1:10){
  1208. curr_data = data.frame(
  1209. score = curr_csv[[metric]][parcel_indices[,i]],
  1210. cohens_d = turb_CohensD[parcel_indices[,i],i],
  1211. lambda = as.factor(rep(lambda_flattened[i],100)),
  1212. parcel = c(1:1000)[parcel_indices[,i]])
  1213. curr_data$group = as.vector(my_parcels[,i])
  1214. data_sensory_assoc = rbind(data_sensory_assoc, curr_data)
  1215. }
  1217. #Grouping by lambda
  1218. data_sensory_assoc = data_sensory_assoc %>%
  1219. group_by(lambda)
  1220. data_sensory_assoc$lambda = as.numeric(as.character(data_sensory_assoc$lambda))
  1222. # Model Fitting
  1224. # 1) Logistic models (modelling the sigmoidal shape)
  1226. # Selection/Computation of initial values for upper, lower asymptode, slope and inflection point
  1227. asymp_up_start = max(data_sensory_assoc$score) # upper asymptode
  1228. asymp_low_start = min(data_sensory_assoc$score) # lower asymptode
  1229. x0_start = mean(range(data_sensory_assoc$lambda)) # inflection point
  1230. k_start = 1 # slope of inflection. negative for decreasing shape
  1232. # Fit 4-parameter logistic model
  1233. logistic_model_4 = nlsLM(
  1234. formula = score ~ asymp_low + (asymp_up - asymp_low) / (1 + exp(k * (lambda - x0))),
  1235. data = data_sensory_assoc,
  1236. start = list(
  1237. asymp_up = asymp_up_start,
  1238. asymp_low = asymp_low_start,
  1239. k = k_start,
  1240. x0 = x0_start
  1241. ),
  1242. trace = TRUE
  1243. )
  1244. print(summary(logistic_model_4))
  1246. #Fitted values for logistic models
  1247. data_sensory_assoc$fitted_log4 = predict(logistic_model_4)
  1249. # 2) Linear Model
  1250. linear_model = lm(score ~ lambda, data = data_sensory_assoc)
  1251. data_sensory_assoc$fitted_linmod = predict(linear_model)
  1253. # 3) Polynomial model
  1254. poly_model = lm(score ~ stats::poly(lambda, 3, raw = TRUE), data = data_sensory_assoc)
  1255. data_sensory_assoc$fitted_poly = predict(poly_model)
  1257. #Comparing the models
  1258. #RSS: raw Model fit., AIC: balance between fit and complexity, BIC: penalizes model complexity more
  1259. c_log_4 = c(sum(residuals(logistic_model_4)^2), AIC(logistic_model_4), BIC(logistic_model_4))
  1260. c_lin = c(sum(residuals(linear_model)^2), AIC(linear_model), BIC(linear_model))
  1261. c_poly = c(sum(residuals(poly_model)^2), AIC(poly_model), BIC(poly_model))
  1262. model_comp = data.frame(
  1263. linear = c_lin,
  1264. logistic_4 = c_log_4,
  1265. polynomial = c_poly
  1266. )
  1267. rownames(model_comp) = c("RSS", "AIC", "BIC")
  1269. gradient_mod_comp[["measures"]] = rownames(model_comp)
  1270. gradient_mod_comp[[names(weighted_sydnor)[metric]]] = model_comp
  1271. gradient_mod_comp[[paste0(names(weighted_sydnor)[metric],"_", "models")]] = list(linear_model,poly_model,logistic_model_4)
  1274. ### Plotting
  1275. # 1 is Control, 0 is TD
  1276. data_sensory_assoc$group = factor(data_sensory_assoc$group,
  1277. levels = c(0, 1),
  1278. labels = c("Autism", "Control"))
  1280. #Compute median for each lambda for plot
  1281. my_medians = data_sensory_assoc %>%
  1282. group_by(lambda) %>%
  1283. summarise(median_score = median(score))
  1285. my_medians$lambda = as.numeric(my_medians$lambda)
  1286. data_sensory_assoc$lambda = as.factor(data_sensory_assoc$lambda)
  1287. data_sensory_assoc$lambda = factor(data_sensory_assoc$lambda, levels = rev(levels(data_sensory_assoc$lambda)))
  1289. #Direction of gradient for that specific metric
  1290. if (gradient_directions$direction[metric]){
  1291. colors = c("goldenrod1","white","#6f1282" )
  1292. } else{
  1293. colors = c("#6f1282","white", "goldenrod1")}
  1295. #Plotting the association of node-turb with the gradient scores for each lambda
  1296. curr_plot = ggplot(aes(x = lambda, y = score), data = data_sensory_assoc) +
  1297. geom_jitter(aes(fill = score, size = cohens_d),
  1298. shape = 21, stroke = 2, colour = "black", width = 0.1, show.legend = c(fill = TRUE, size = FALSE)) +
  1299. labs(size = "Cohen's d") +
  1300. scale_fill_gradientn(colors = colors, guide = guide_colorbar(barwidth = 1.5, barheight = 10)) +
  1301. scale_size_continuous(range = c(-4,10)) +
  1302. geom_boxplot(aes(x = lambda, y = score), size = 1.5, colour = "black", fill = NA,
  1303. position = position_dodge(width = 1), outlier.shape = NA) +
  1304. ylab(names(weighted_sydnor)[metric]) +
  1305. xlab("Lambda (in 1/mm)") +
  1306. theme_minimal() +
  1307. theme(
  1308. axis.title = element_blank(),
  1309. axis.text = element_text(size = 30),
  1310. plot.title = element_blank(),
  1311. axis.title.x = element_blank(),
  1312. axis.text.x = element_blank(),
  1313. axis.ticks.x = element_blank(),
  1314. legend.title = element_blank(),
  1315. legend.text = element_text(size = 30)
  1316. ) +
  1317. geom_line(data = data_sensory_assoc, aes(x = lambda, y = fitted_log4, group = 1), color = "black", linetype = "solid", linewidth = 4)
  1318. print(curr_plot)
  1319. #ggsave(paste0(figures_directory,"/Figure_5/",names(weighted_sydnor)[metric],"_S_SA.jpg"), curr_plot, width = 13.4, height = 8.4, dpi = 300)
  1320. }
  1321. ```
  1323. 3.3 Printing the association of Node-Turbulence to Sydnor gradients
  1325. ```{R}
  1326. #Loading data. The computation of the effect sizes of the 400 parcels has been done with the above python script and is only imported here. The matrices have the size (400x10), where each row corresponds to a parcel and each column to a Lambda value. The matrices are filled with Cohen's d effect sizes (TD-autism)
  1328. SA_data_path = paste0(root,"/data/gradient_analysis/")
  1330. #Loading the converted parcel-level data (for Schaefer400).
  1331. parcel_names = c("parcel_turb_autism_s400", "parcel_turb_TD_s400", "parcel_turb_cohensD_s400")
  1332. parcel_data = list()
  1334. for (currName in parcel_names){
  1335. parcel_data[[currName]] = data.frame(read.csv(paste0(SA_data_path,currName,".csv")))
  1336. }
  1338. # Load schafer atlas
  1339. data("schaefer17_400")
  1340. df = schaefer17_400$data
  1341. my_labs = c(schaefer17_400_3d$ggseg_3d[[1]]$label[-1], schaefer17_400_3d$ggseg_3d[[2]]$label[-1])
  1342. sens_assoc = read_csv(paste0(root,"/code/Sensorimotor_Association_Axis_AverageRanks.csv"))
  1343. brainmaps_gradients = read_csv(paste0(root,"/code","/brainmaps_schaefer.csv"))
  1345. #Preallocation
  1346. curr_lambda_plots = list()
  1347. for (currLam in 1:10){
  1348. #plot init
  1349. curr_data = data.frame(
  1350. dat = parcel_data$parcel_turb_cohensD_s400[,currLam])
  1351. df_gradients = cbind(curr_data, brainmaps_gradients, sens_assoc)
  1353. # Step 1: Filter atlas to labels in my_labs
  1354. df_filtered = df %>%
  1355. dplyr::filter(label %in% my_labs)
  1357. # Step 2: Match each label with its corresponding value
  1358. df_ordered = df_filtered %>%
  1359. left_join(df_gradients, by = "label")
  1361. # Step 3: rearrange the atlas
  1362. schaefer17_400$data = df_ordered
  1364. colors = c("goldenrod1","white","#6f1282" )
  1366. #Adjust color scale
  1368. # Map to quantiles for smoother gradient
  1369. probs = seq(0, 1, length.out = 100)
  1370. quant_vals = quantile(df_ordered$dat, probs, na.rm = TRUE)
  1372. # curr_plot
  1373. curr_plot = ggseg(
  1374. atlas = schaefer17_400,
  1375. mapping = aes(fill = dat),
  1376. color = "black",
  1377. size = 0.3
  1378. ) +
  1379. scale_fill_gradientn(
  1380. colors = colors,
  1381. values = scales::rescale(quant_vals, from = range(quant_vals)),
  1382. limits = range(quant_vals),
  1383. oob = scales::squish
  1384. ) +
  1385. theme_void() +
  1386. theme(
  1387. legend.title = element_blank(),
  1388. legend.text = element_blank(),
  1389. legend.key.height = unit(0.4, "cm"),
  1390. legend.key.width = unit(0.4, "cm"),
  1391. legend.position = "none"
  1392. )
  1393. #Add to list
  1394. curr_lambda_plots[[currLam]] = curr_plot
  1395. }
  1397. # S–A gradient
  1398. SA_plot = ggseg(
  1399. atlas = schaefer17_400,
  1400. mapping = aes(fill = finalrank.wholebrain),
  1401. color = "black",
  1402. size = 0.3) +
  1403. scale_fill_gradientn(colors = colors, values = c(0, 0.5, 1)) +
  1404. theme_void() +
  1405. theme(legend.position = "none")
  1407. #Functional gradient
  1408. FUNC_plot = ggseg(
  1409. atlas = schaefer17_400,
  1410. mapping = aes(fill = G1.fMRI),
  1411. color = "black",
  1412. size = 0.3) +
  1413. scale_fill_gradientn(colors = colors, values = c(0, 0.5, 1)) +
  1414. theme_void() +
  1415. theme(legend.position = "none")
  1417. #Plotting
  1418. plot_selected_lambdas = curr_lambda_plots[[1]]/curr_lambda_plots[[4]]/curr_lambda_plots[[7]]/curr_lambda_plots[[10]]/SA_plot/FUNC_plot
  1419. print(plot_selected_lambdas)
  1420. #ggsave(paste0(figures_directory,"/Figure_5/parcel_turb_S400.jpg"), plot_selected_lambdas, width = 16, height = 16, dpi = 300)
  1422. #Heatmaps of correlation values for SA and func gradient with parcel-level turbulence
  1423. #load rho values and p-values
  1424. corr_results = read.csv(paste0(root, "/data/gradient_analysis/corr_results.csv"))
  1426. #DF for plotting
  1427. df_SA_heatmap = data.frame(
  1428. Variable = rep(factor(Lambda), 2),
  1429. rho = round(c(corr_results$rho_SA, corr_results$rho_func),2),
  1430. neglog10p = c(-log10(corr_results$p_SA), -log10(corr_results$p_func)),
  1431. x = rep(c(0,0.6), each = 10)
  1432. )
  1434. #Create the colormap
  1435. SA_heatmaps = ggplot(df_SA_heatmap, aes(x = x, y = Variable, fill = neglog10p)) +
  1436. geom_tile(color = "white", width = 0.6) +
  1437. geom_text(aes(label = sprintf("%.2f", rho)), size = 10) +
  1438. scale_fill_distiller(palette = "OrRd", direction = 1) +
  1439. scale_x_discrete(expand = c(0, 0.1)) +
  1440. labs(
  1441. x = "",
  1442. y = "",
  1443. fill = "-log10(p)"
  1444. ) +
  1445. theme_minimal() +
  1446. theme(
  1447. axis.text.x = element_blank(),
  1448. axis.ticks.x = element_blank(),
  1449. axis.text.y = element_blank(),
  1450. axis.ticks.y = element_blank(),
  1451. panel.grid = element_blank(),
  1452. legend.position = "none",
  1453. plot.margin = margin(5, 5, 5, 5)
  1454. ) +
  1455. coord_fixed(ratio = 0.25)
  1457. print(SA_heatmaps)
  1458. ##ggsave(paste0(figures_directory,"/Figure_5/SA_heatmaps.jpg"), SA_heatmaps, dpi = 300)
  1459. ```

Turbulence_postproc_analysis.Rmd at commit 15de1d4, no license · at the source

Overview

Authors: Benedikt P Fuhr1,2, Yonatan Sanz Perl3,4, Ines Severino1,2, Morten L Kringelbach5,6,7,8, Gustavo Deco3,4,6, Henricus G Ruhé9,10, Michael V Lombardo1
ORCID iDs: Benedikt P Fuhr, Michael V Lombardo
  1. Laboratory for Autism and Neurodevelopmental Disorders, Center for Neuroscience and Cognitive Systems, Istituto Italiano di Tecnologia, Rovereto, Italy
  2. Center for Mind/Brain Sciences, University of Trento, Rovereto, Italy
  3. Center for Brain and Cognition, Computational Neuroscience Group, Faculty of Medicine and Life Sciences, Universitat Pompeu Fabra, Barcelona, Spain
  4. Institució Catalana de la Recerca i Estudis Avançats (ICREA), Barcelona, Spain
  5. Centre for Eudaimonia and Human Flourishing, Linacre College, University of Oxford, Oxford, United Kingdom
  6. International Centre for Flourishing, Universities of Oxford (United Kingdom), Aarhus (Denmark), and Pompeu Fabra (Spain)
  7. Department of Clinical Medicine, Aarhus University, Aarhus, Denmark
  8. Department of Psychiatry, University of Oxford, Oxford, United Kingdom
  9. Department of Psychiatry, Radboudumc, Nijmegen, the Netherlands
  10. Donders Institute for Brain, Cognition and Behavior, Radboud University, Nijmegen, the Netherlands
Institutions: University of Trento (Italy); Italian Institute of Technology (Italy); Center for Neuroscience and Cognitive Systems (Italy); Institució Catalana de Recerca i Estudis Avançats (Spain); Universitat Pompeu Fabra (Spain); Aarhus University (Denmark); University of Oxford (United Kingdom); Radboud University Nijmegen (Netherlands); Radboud University Medical Center (Netherlands); Donders Institute for Brain, Cognition and Behaviour (Netherlands)
Journal: Biological psychiatry global open science, volume 6, issue 5, article 100787
Dates: received 20 March 2026; accepted 26 June 2026; published online 8 July 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.bpsgos.2026.100787 · PMID 42602560 · PMCID PMC13475223 · OpenAlex W4417290912
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.64898/2025.12.11.25342053
Status: code verified
Categories: fMRI (modality), autism (population), systems (subfield)
Methods: Statistics, Spectral & time-frequency, fMRI & imaging
Keywords: Autism, Brain dynamics, Cortical gradients, Functional connectivity, Imaging, Resting-state fMRI
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: European Research Council
Citations: not cited yet (Europe PMC); 65 references in the paper

Abstract

Background: Prevailing theories propose that autism is characterized by local cortical overconnectivity and long-range underconnectivity, but empirical evidence remains mixed.

Methods: Here, we applied the turbulence dynamics framework to the ABIDE (Autism Brain Imaging Data Exchange) dataset (N = 1009) to examine how functional synchronization profiles dynamically change over time and across the cortex over different spatial scales.

Results: Autistic individuals showed increased short-range and reduced long-range functional synchronization variability over time, as well as reduced synchronization strength across all spatial scales. Synchronization also decayed more rapidly with distance and exerted weaker influence across scales in autism. These distance-specific alterations suggest that local hyperconnectivity may be associated with turbulent synchronization dynamics that fail to propagate coherently across the cortex, resulting in an overly rigid brain organization at longer distances. Mapping these effects onto the sensorimotor–association cortical gradient revealed increased variability in sensorimotor regions and decreased variability in the association cortex.

Conclusions: Together, we found evidence of disturbances in functional synchronization dynamics at different spatial scales and along hierarchical brain gradients in autistic individuals. These results consolidate ideas about dynamic functional connectomic organization in autism and situate these alterations along hierarchical brain gradients that are closely linked to neurodevelopmental processes.

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

Repository

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

gitlab.iit.it/bmp006-public/land_iit/turbulence_abide

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 15de1d443f70190d86a707b7ff8543f64b2f26b3, 21 March 2026
Languages: Python (6), MATLAB (3), Shell (3), Jupyter (2), R (1)
Size: 2,032 files, 15 scripts
Software Heritage: not archived
Found in: the acknowledgements
Holds: README, environment (requirements.txt), 2 notebooks
Not found: license file, CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (5 files), pandas (5 files), AFNI (3 files), BrainSpace (3 files), SciPy (3 files), Statistics and Machine Learning Toolbox (2 files), Matplotlib (2 files), FSL (1 file), lmerTest (1 file), Optimization Toolbox (1 file), Signal Processing Toolbox (1 file), NiBabel (1 file), tidyverse (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
16 files
At the source: gitlab.iit.it/bmp006-public/LAND_IIT/turbulence_abide/

Tracing map

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

What the map holds:

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

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

Data

No dataset and no data link were found in the paper.

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 2, 28 September 2026

  • Authors: added Benedikt P Fuhr (0009-0005-5879-9198); Michael V Lombardo (0000-0001-6780-8619); removed Benedikt P Fuhr; Michael V Lombardo

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 6 keywords, 1 funder, 62 references.

Cite

This paper

Fuhr, B. P., Perl, Y. S., Severino, I., Kringelbach, M. L., Deco, G., Ruhé, H. G., & Lombardo, M. V. (2026). Dynamic Functional Synchronization Profiles in Autism Differ by Spatial Scale and Along Hierarchical Cortical Gradients. Biological psychiatry global open science, 6(5), 100787. https://doi.org/10.1016/j.bpsgos.2026.100787

BibTeX

@article{fuhr2026dynamic,
author = {Fuhr, Benedikt P and Perl, Yonatan Sanz and Severino, Ines and Kringelbach, Morten L and Deco, Gustavo and Ruhé, Henricus G and Lombardo, Michael V},
title = {{Dynamic Functional Synchronization Profiles in Autism Differ by Spatial Scale and Along Hierarchical Cortical Gradients}},
journal = {Biological psychiatry global open science},
year = {2026},
month = jul,
volume = {6},
number = {5},
pages = {100787},
publisher = {Elsevier},
issn = {2667-1743},
doi = {10.1016/j.bpsgos.2026.100787},
url = {https://doi.org/10.1016/j.bpsgos.2026.100787},
pmid = {42602560},
pmcid = {PMC13475223}
}

RIS

TY - JOUR
AU - Fuhr, Benedikt P
AU - Perl, Yonatan Sanz
AU - Severino, Ines
AU - Kringelbach, Morten L
AU - Deco, Gustavo
AU - Ruhé, Henricus G
AU - Lombardo, Michael V
TI - Dynamic Functional Synchronization Profiles in Autism Differ by Spatial Scale and Along Hierarchical Cortical Gradients
T2 - Biological psychiatry global open science
J2 - Biol Psychiatry Glob Open Sci
PY - 2026
DA - 2026/07/08
VL - 6
IS - 5
SP - 100787
SN - 2667-1743
PB - Elsevier
DO - 10.1016/j.bpsgos.2026.100787
UR - https://doi.org/10.1016/j.bpsgos.2026.100787
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.bpsgos.2026.100787",
"type": "article-journal",
"title": "Dynamic Functional Synchronization Profiles in Autism Differ by Spatial Scale and Along Hierarchical Cortical Gradients",
"container-title": "Biological psychiatry global open science",
"author": [
{
"family": "Fuhr",
"given": "Benedikt P"
},
{
"family": "Perl",
"given": "Yonatan Sanz"
},
{
"family": "Severino",
"given": "Ines"
},
{
"family": "Kringelbach",
"given": "Morten L"
},
{
"family": "Deco",
"given": "Gustavo"
},
{
"family": "Ruhé",
"given": "Henricus G"
},
{
"family": "Lombardo",
"given": "Michael V"
}
],
"container-title-short": "Biol Psychiatry Glob Open Sci",
"volume": "6",
"issue": "5",
"page": "100787",
"DOI": "10.1016/j.bpsgos.2026.100787",
"PMID": "42602560",
"PMCID": "PMC13475223",
"ISSN": "2667-1743",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.bpsgos.2026.100787",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
8
]
]
}
}

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

Similar papers

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

[1] doi:10.1038/s41593-026-02287-z [code]
Autism subtypes identified using cross-species functional connectivity analyses.
Journal: Nature neuroscience
In common: AFNI, FSL, NiBabel, 2 other tools, autism, fMRI, 12 references, author Michael V Lombardo
[2] doi:10.1038/s41467-026-76011-7 [code]
Human cortex organizes dynamic co-fluctuations along the sensorimotor-association axis.
Journal: Nature communications
In common: BrainSpace, AFNI, FSL, 7 other tools, systems, 6 references
[3] doi:10.1038/s42003-026-10276-y [code]
The cellular correlates and adolescent reorganisation of cortical myelination networks in the common marmoset.
Journal: Communications biology
In common: BrainSpace, AFNI, Optimization Toolbox, 8 other tools, 4 references
[4] doi:10.1038/s41531-026-01354-3 [code]
Neuromodulation-induced normalization of cortical metastable dynamics signatures in Parkinson's disease.
Journal: NPJ Parkinson's disease
In common: BrainSpace, Signal Processing Toolbox, NiBabel, 5 other tools, 8 references
[5] doi:10.1038/s41467-026-73668-y [code]
Convergent and divergent brain-cognition development in early adolescence.
Journal: Nature communications
In common: AFNI, FSL, Signal Processing Toolbox, 7 other tools, fMRI, 6 references
[6] doi:10.1038/s41467-026-71270-w [code]
Spatiotemporal dynamics of the human cortical functional hierarchy across the lifespan.
Journal: Nature communications
In common: BrainSpace, FSL, NiBabel, 5 other tools, fMRI, 6 references
[7] doi:10.64898/2026.03.12.710517 [code]
Cortical excitability inversely modulates fMRI connectivity via low-frequency neuronal coupling
Journal: bioRxiv (preprint)
In common: AFNI, Optimization Toolbox, FSL, 6 other tools, fMRI, systems, 3 references
[8] doi:10.1038/s41586-026-10631-3 [code]
A prognostic human brain network for diffuse midline glioma.
Journal: Nature
In common: BrainSpace, Optimization Toolbox, FSL, 7 other tools, 3 references
[9] doi:10.1038/s41593-026-02205-3 [code]
Competitive interactions shape mammalian brain network dynamics and computation.
Journal: Nature neuroscience
In common: FSL, Signal Processing Toolbox, NiBabel, 5 other tools, 5 references
[10] doi:10.1038/s41467-026-75959-w [code]
Charting higher-order models of brain function beyond pairwise interactions.
Journal: Nature communications
In common: BrainSpace, NiBabel, Statistics and Machine Learning Toolbox, 4 other tools, 7 references

Contribute

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

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

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

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