4.1 Participants and procedures
Data were taken from NeuroIMAGE II, the third wave of an integrated genetics-cognition-MRI-phenotype project on ADHD 38. It includes resting-state fMRI data of 249 individuals with ADHD as well as age- and sex-matched healthy controls. Initial recruitment criteria for ADHD participants were an ADHD combined type diagnosis, availability of one or more siblings, age between 6–18 years and availability of subject, sibling, and at least one biological parent for DNA collection. Exclusion criteria (for all participants) were IQ (as measured by Wechsler Intelligence Scale for Children/Wechsler Adult Intelligence Scale) < 70, diagnoses of autism or schizophrenia, and neurological disorders. For controls, it was required that neither they nor any of their first-degree relatives had a prior ADHD diagnosis.
Reassessment of the ADHD diagnosis was established by combining information from the Kiddie Schedule for Affective Disorders (K-SADS) 39 and parent, teacher, and self-report versions of the Conners’ rating scale (CPRS-R:L, CTRS-R:L & CAARS-R:L) 40–42. All diagnostic and phenotypic data was acquired on the same day as the fMRI data. A detailed description of the diagnostic procedures is given by von Rhein et al. 38.
Participants with a subthreshold ADHD diagnosis (2–4 ADHD-specific symptoms), left-handedness, excessive movement during the scanning, or insufficient quality of rs-fMRI or questionnaire data were excluded. 147 participants were used for the final analysis, including 31 participants with ADHD but without hyperactivity-impulsivity symptoms, i.e., predominantly inattentive ADHD (ADHD-I), 25 participants with ADHD and hyperactivity-impulsivity symptoms (ADHD-C/H), i.e., 21 with combined type ADHD and 4 with predominantly hyperactive-impulsive ADHD, and 91 healthy controls. Demographic information are given in Table 1.
Forty-eight hours prior to testing, stimulant medication use was discontinued. Data acquisition took place at the Donders Institute for Cognitive Neuroimaging, Radboud University Nijmegen, Netherlands. Participants (and their parents when < 18 years old) gave written informed consent for participation. Ethical approval was granted by the regional ethics board (Centrale Commissie Mensgebonden Onderzoek: CMO Regio Arnhem Nijmegen, ABR: NL41950.091.12). Data analysis was pre-registered using the open science framework (osf.io/rdyp6) 43.
4.2 Resting-state fMRI data acquisition and preprocessing
Whole-brain imaging was performed on a 1.5 T Magnetom Avanto (Siemens AG, Erlangen, Germany). BOLD-sensitive resting-state functional volumes were acquired using a T2*-weighted EPI sequence (TR = 1960 ms, TE = 40 ms). Each of the 266 volumes consisted of 37 axial slices of size 64x64 (flip angle = 80°, FoV = 224 x 224 mm2, voxel- size = 3.5 x 3.5 x 3.0 mm3, inter-slice gap = 0.5 mm). T1-weighted high-resolution structural volumes were acquired with an MPRAGE sequence (TR = 2730 ms, TE = 2.95 ms, TI = 900 ms, flip angle = 9°, FoV = 256 x 256 mm2, voxel- size = 1.0 x 1.0 x 1,0 mm3, GRAPPA 2).
Preprocessing mostly relied on FMRIB algorithms (FSL 5.0.11; Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012). The resting-state time series data were skull stripped, realigned to the middle volume of the series, co-registered to the structural T1, and spatially smoothed using a 6 mm full width at half maximum Gaussian kernel (FWHM). ICA-AROMA 45 was used to account for secondary movement artefacts. Residual noise was further reduced by nuisance regression including a linear trend and average times series measured within the white matter and cerebrospinal fluid. High-pass filtering was conducted at 0.008 Hz. Prior to network analysis, time series were warped to MNI152 space (Montreal Neurological Institute, Montreal, Canada). Root mean squared framewise displacement was calculated. A threshold of 0.25 was applied to exclude 25 participants with extreme movement from further analysis. Root mean squared framewise displacement did not significantly differ between the groups (healthy controls: 0.087 ± 0.071; ADHD-I: 0.128 ± 0.096; ADHD-C/H: 0.111 ± 0.096).
4.3 Graph analysis
For graph analysis we used Python 3.5 with NetworkX 46. Parcellation of preprocessed time series data was realized using a hemisphere-specific functional brain template with 268 parcels 47. It was created using graph-theoretically based parcellation that ensures functional homogeneity within the parcels of the atlas, even across different individuals. The atlas thus minimizes the likelihood that different functional areas lie within a single parcel 47. For 221 relevant parcels, covered by the MR measurement, subject-specific average intensity time series were calculated. Correlation matrices were created by computing pairwise Pearson’s correlations between the extracted time series. Matrices were Fisher’s z-transformed and absolute values were taken. Absolute value transformation was performed as preprocessing of the present data did not involve global signal regression and as anti-correlations are thought to be functionally relevant 48. Due to the naturally low density of negative correlations, we refrained from performing specific analysis for positive and negative correlations. To distinguish differences in network density from those of network topology 49, matrices were binarized based on seven equally spaced density thresholds with a minimum density of 0.10 and maximum density of 0.40 50. In this range of low to medium network densities, previous studies found significant associations between network topology and ADHD symptoms 11,12,14. Thus, seven threshold-specific graphs for each subject were investigated in the following graph analysis. Nodal topology measures were calculated for 46 of the 221 nodes. These 46 nodes were chosen based on their association with parcels overlapping with the orbitofrontal cortex, ventromedial prefrontal cortex, anterior cingulate cortex, insula, ventral striatum, amygdala, and hippocampus as defined by the Harvard-Oxford Brain Atlas by more than 30%. These brain regions have been previously documented to be involved in emotion processing, its implicit regulation, and brain dysfunctions in ADHD 7. Since we focused on implicit regulation, regions more strongly related to explicit processes, such as cognitive reappraisal, (e.g., frontoparietal network) were not considered.
Here, we captured the centrality of each node and analyzed whether nodes were highly integrated and connected, either locally towards their direct neighboring nodes or globally with the entire network. Thus, our focus is on six nodal measures, that is betweenness, closeness, eigenvalue centrality, clustering coefficient, nodal efficiency, and local efficiency, which were used in previous studies on ADHD and showed ADHD-specific deviations 11,12,14. Betweenness, closeness, and eigen-value centrality are measures that describe the centrality of a node within a network. While betweenness describes how often a node is part of the shortest path between two other nodes, closeness describes how many of the theoretically possible direct connections to other nodes actually exist. Eigenvalue centrality also considers the centrality of the node’s direct neighbors. The clustering coefficient of a node describes how strongly its neighboring nodes are interconnected. Efficiency values indicate how directly nodes can be reached from other nodes of the network. In the case of nodal efficiency, this refers to the shortest connection from a particular node to all other nodes, while in the case of local efficiency it refers to the efficiency amongst the nodes adjacent to a particular node of interest. All measures entered into separate SEM models described in 2.4.
Density-integrated topology measures were calculated 49. Differences between populations of weighted networks may be due to the networks’ wiring costs and not the targeted topological features. Density-integration of measures from binarized networks, however, can eliminate cost-related differences and also allow the assessment of topology measures under different density-thresholds. Figure 2 summarizes the functional connectivity and network analysis.
4.4 Statistical analysis and structural equation modeling
Variable Selection. Statistical analyses were conducted with R software 51. Lavaan was used for SEM 52. Prior to SEM, multi-group confirmatory factor analysis was used to identify data suitable for calculation of the latent emotion dysregulation variable. NeuroIMAGE II includes 6 tasks, questionnaires, or questionnaire subscales, respectively, that gauge emotional problems, emotional lability, and associated features: the Kessler Psychological distress scale, NESDA version of K-10 53 with 10 items for the assessment of emotional problems including anxiety and depression, the strength and difficulties questionnaire subscales (SDQ; five items about anxieties, worries, happiness, and physical symptoms of emotional stress for the emotional symptoms subscale & five items about temper tantrums, compliance, quarrelsomeness, stealing, and lying for the conduct symptoms subscales) 54, the emotional lability subscale of the Conners’ parent rating scale (CPRS-R:L consisting of three items for unpredictable mood changes, temper tantrums, and tearfulness) 42, the Inventory of Callous-Unemotional traits (ICU) 55 (24 items with a callousness and unemotional traits score), and the MINDS Testmanager’s gradual emotion recognition task (GERT) 56 (accuracy of correct emotion classification).
Mediation structural equation model. The latent variable was used to investigate the relationship between functional brain network activity, emotion dysregulation, and ADHD. In the initial, preregistered analysis plan, we intended to use SEM mediation models in which the relationship between topological measures and the emotional latent variable was mediated by ADHD scores (CPRS-R:L) without taking into account the different ADHD subtypes. However, those models did not produce good model-fit (as measured by the standard goodness-of-fit measures described below) and had no significant results. One reason may be that the ADHD scores were derived from a parents’ questionnaire for children’s ADHD symptoms (CPRS-R:L) which may be less valid than a clinical diagnosis.
Three-group structural equation model. Alternatively, we chose a multi-group SEM approach using the clinical diagnoses. The diagnoses reflect all available diagnostic information and may thus provide more accurate models. We investigated the group dependent differences (i.e., healthy controls, ADHD-I participants, & ADHD-C/H participants) in the association between topology measures and the latent variable gauging emotion dysregulation. For each node of interest and each density-specific topology measure, one multi-group model was built. Models consisted of the latent variable, questionnaire scores, associated parameter estimates and variances as well as the nodal topology measure variable with its regression parameter for the latent variable (see Fig. 3). Factor loadings of the latent variable were fixed across groups, latent variable variance was standardized, and to account for between-group mean differences, group-specific intercepts were added. Models were compared with almost identical models, in which however the regression parameters between the network topology variable and the latent emotion dysregulation variable were fixed across groups. Model estimation was conducted using maximum likelihood procedures. To evaluate the significance of the group effects, model-specific χ2-values were compared between the two models (χ2-difference-test). Following the approach of Ginestet et al. 49, a model was considered significant if it revealed a significant effect on the density-integrated level (averaged over densities). Subsequently, significance was tested on each density threshold to provide detailed information on whether effects were found with stronger or weaker network connections.
Post-hoc analysis, multiple comparison procedures and goodness-of-fit. We used post-hoc two-group SEMs with Bonferroni-corrected p-values to investigate pairwise between-group differences. These two-group models resembled the three-group models described above, except that each only considered participants from two of the three groups. Alpha inflation due to multiple comparisons was controlled by the Benjamini-Hochberg false discovery rate (FDR) procedure 57. The comparative fit index (CFI), standardized root mean square residual (SRMR), and root mean square error of approximation (RMSEA) were calculated to evaluate the models’ goodness-of-fit. Acceptable goodness-of-fit measures should at least be above 0.95 for the CFI, below 0.08 for the SRMR, and below 0.10 for the RMSEA 58.
Additional analyses. To further control the robustness of the results, we evaluated whether group-specific differences in overall connectivity strength of the underlying adjacency matrices exist 59. Furthermore, an alternative parcellation scheme was used to investigate the robustness of the results with respect to the chosen parcellation approach. For this, we used the AAL atlas 60, further subdivided with K-mean clustering to adapt the size and number of parcels to that of Finn et al. (2015). Procedures for calculating topology measures and SEM analyses were conducted as described above.
Data availability statement
The data are property of the Donders Institute for Cognitive Neuroimaging and are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions. Further information about the NeuroIMAGE project may be requested via the project leader Jan K. Buitelaar.