Analysis of MRI data

MRI data were converted to the Brain Imaging Data Structure (BIDS) [47] using in-house scripts preceding the BiDirect-BIDS-ConverteR [48, 49]. Facial features were removed from the T1-weighted anatomical data [50]. Main MRI data pre-processing was carried out using fMRIPrep [51] (version 20.0.7, RRID:SCR_016216) briefly consisting of motion estimation and correction, co-registration of fMRI and structural MRI data, estimation of noise regressors, and standard space normalization. Please consult the Supplementary Methods for further details. Subsequent actual denoising was carried out using fMRIDenoise [52] (version 0.2.1), comprising regressing out 24 head motion parameters (3 translations, 3 rotations, their 6 temporal derivatives, and their 12 quadratic terms) [53], 8 physiological noise parameters (mean physiological signals from white matter and cerebrospinal fluid, their 2 temporal derivatives, and 4 quadratic terms) [53] as well as movement spike regression based on frame-wise displacement (FD > 0.5 mm) and so-called “DVARS” (> 3) thresholds [54], temporally filtering (0.008–0.08 Hz), and, finally, smoothing the resulting standard space image with a Gaussian kernel (FWHM = 6 mm). This pre-processing leads to denoised fMRI data in a common standard space as input for further analyses.

The following steps were taken for MRI data quality control (numbers of excluded participants reported in the section “participants”): Structural MRI data were screened by a radiologist for incidental findings and major artifacts. fMRIPrep reports were reviewed for registration errors and image artifacts (including signs of strong motion artifacts in the carpet plots). Subject exclusion for excessive head motion (see section “Participants”) was based on pre-processing criteria (mean frame-wise displacement, FD > 0.3 mm or maximum FD > 5 mm or more than 20% outlier data points). FD did not differ significantly between FAS and controls. However, there was a trend towards higher mean and maximum FD in the FAS group (see Supplementary Table 1).

Functional connectivity analyses were based on a cortical atlas (“Schaefer atlas”) derived from rs-fMRI data in 1489 participants. The atlas was obtained from TemplateFlow (RRID:SCR_021876) to match the dimensions of the fmriprep outputs [55]. The atlas version with 400 parcels adopted here is the most extensively validated version of this atlas, e.g., regarding stability and correspondence with markers of brain function [56]. The individual parcels in the published atlas have been matched to 17 non-overlapping networks from a previously established atlas by Yeo et al. [57]. Ten cognition-related components out of these 17 networks were selected for further analysis: the dorsal attention network (2 sub-networks A and B), the salience / ventral attention network (2 sub-networks A and B), the mainly fronto-parietal control network (3 sub-networks A-C), and the default mode network (3 sub-networks A-C). Time-series extraction from the pre-processed fMRI data and calculation of z-transformed Pearson correlation coefficients as primary measures of FC were carried out with the Data Processing Assistant for Resting-State fMRI (DPARSF, version 5.2, RRID:SCR_002372) [58] based on Matlab 2019b (The MathWorks, Natick, MA, USA, RRID:SCR_001622). As a basis for subsequent modelling, we carried out multiple linear models (one model per pair of regions), comparing z-transformed correlation coefficients (dependent variable) among regions of interest between FAS participants and controls. The following independent variables were included in the multiple linear regression models: group (FAS vs. control, categorical), sex (categorical), age (normalized to center: 0, and standard deviation: 1), mean FD, and a constant term (intercept). Main effects of the factor group (FAS vs. controls, two-sided) are the basis of the subsequent analyses at different spatial scales. A previous analogous analysis based on multiple two-tailed two-sample t-tests, i.e., not including age, sex, and head movement with similar results has been made available in the preprint version of this article [59].

We followed a hierarchical statistical testing approach with three levels of analysis: (1) first, determining whether there is any alteration of FC in FAS participants compared with controls in the full connectome of 243 regions constituting these 10 cognition-related networks (i.e., an omnibus test) globally and subsequently aiming to resolve these findings, (2) to the level of either FC within each network or between-network connectivity, and finally (3) to individual connections.

The omnibus test on the full connectome (main effect of group: FAS vs. control) is based on the “Higher Criticism” (HC) approach [60] in an improved version [61] implemented in Matlab (https://www.stat.cmu.edu/~jiashun/Research/software/HC/). HC statistics can be applied in order to test whether there are any non-zero effects within a large number of individual tests carried out in high-dimensional datasets. They are thus suitable to identify the existence of rare and weak effects in such data [60]. HC follows the rationale of p-value histogram analyses: Under a null-hypothesis of only zero effects in multiple parallel tests, an equal distribution of p-values is expected. Under the alternative hypothesis of existing non-zero effects, there is an excess of low p-values [62]. In simplified terms, HC statistics quantitatively test a joint hypothesis of such an excess of low p-values [60, 62]. HC has been increasingly popular for detecting effects in high-dimensional data such as in genetic [60] and economic [63] research. It has been argued that HC could be favorable for the detection of rare events compared with conventional false discovery rate (FDR) or family-wise error rate (FWE) correction methods [60]. Considering the similarly high dimensionality of FC datasets, HC has recently been applied to rs-fMRI analyses [64, 65]. Both, in this global analysis and the subsequent network-wise HC-based hypothesis tests, we used p-value histograms (main effect of group) as a plausibility check.

Subsequently, we aimed to determine which of the 10 cognition-related networks were affected by within-network functional connectivity differences between FAS and control participants. We therefore carried out equivalent HC tests separately for these networks. Each set of tests included the full set of correlation coefficients between all regions within each network.

For a similar analysis of between-network FC (i.e., to determine whether any between-network FC differences were present), we concatenated all parcels for each of the 10 networks separately. This resulted in a single mask for each network before time-series data extraction. The resulting correlation coefficients were assessed with an equivalent HC test.

In a third level, we aimed to identify single between- and within-network connections exhibiting statistically significant FC differences between FAS participants and controls. Therefore, in contrast to the previously described HC-based joint hypothesis tests, we now FDR-adjusted [66] the individual hypothesis tests (main effect of group: FAS vs. control) of between- and within-network connectivity (q < 0.05), using an FDR implementation in Matlab (https://brainder.org/2011/09/05/fdr-corrected-fdr-adjusted-p-values). This was carried out separately for either all 45 between-network connections or all individual connections, and also within each separate network.

Beyond the static FC analysis, we carried out an exploratory dynamic FC analysis, using a sliding window approach with the DynamicBC toolbox (version 2.2) [67]. FC between all 243 regions in the cognition-related networks was calculated for each individual subject by Pearson linear correlation separately in overlapping windows with a length of 18 consecutive functional volumes equivalent to 45 s, similar to window lengths in previous studies [68, 69], and with a 60% overlap. The resulting time-resolved FC estimates from individual time windows were grouped by similarity (K-means cluster analysis, distance measure: correlation) in order to derive presumed FC states in the entire sample. Established methods were used to estimate the optimal number of clusters and assess the goodness of fit of the clustering solutions: (1) The optimal number of clusters (search range: 2 to 10) was estimated using the Calinski-Harabasz [70] and Davies–Bouldin [71] indices resulting in 2 clusters (see Supplementary Fig. 1). (2) Cluster-separability was estimated by a silhouette analysis. Briefly, the silhouette analysis assesses, how similar an individual element is to other elements within its own cluster compared with elements in other clusters [72]. The following summary measures describing the temporal dynamics of putative FC states were calculated for individual participants: number of transitions (NT) between connectivity states, mean dwell time (MDT) per cluster, and frequency of observing each cluster (FRC). Participants with FAS and controls were compared regarding NT and MDT using t-tests or Mann–Whitney U tests.