Widespread effects of dMRI data quality on diffusion measures in children.

Diffusion magnetic resonance imaging (dMRI) datasets are susceptible to several confounding factors related to data quality, which is especially true in studies involving young children. With the recent trend of large-scale multicenter studies, it is more critical to be aware of the varied impacts of data quality on measures of interest. Here, we investigated data quality and its effect on different diffusion measures using a multicenter dataset. dMRI data were obtained from 691 participants (5-17 years of age) from six different centers. Six data quality metrics-contrast to noise ratio, outlier slices, and motion (absolute, relative, translation, and rotational)-and four diffusion measures-fractional anisotropy, mean diffusivity, tract density, and length-were computed for each of 36 major fiber tracts for all participants. The results indicated that four out of six data quality metrics (all except absolute and translation motion) differed significantly between centers. Associations between these data quality metrics and the diffusion measures differed significantly across the tracts and centers. Moreover, these effects remained significant after applying recently proposed harmonization algorithms that purport to remove unwanted between-site variation in diffusion data. These results demonstrate the widespread impact of dMRI data quality on diffusion measures. These tracts and measures have been routinely associated with individual differences as well as group-wide differences between neurotypical populations and individuals with neurological or developmental disorders. Accordingly, for analyses of individual differences or group effects (particularly in multisite dataset), we encourage the inclusion of data quality metrics in dMRI analysis.

INTRODUCTION

Diffusion magnetic resonance imaging (dMRI) is a noninvasive neuroimaging tool for probing the microstructural architecture of the brain (Basser, Mattiello, & LeBihan, 1994; Merboldt, Hanicke, & Frahm, 1985). In recent years, several dMRI‐based models like Diffusion Tensor Imaging (DTI; Basser, Mattiello, & LeBihan, 1994; Pierpaoli, Jezzard, Basser, Barnett, & Di Chiro, 1996), Diffusion Kurtosis Imaging (DKI; Jensen & Helpern, 2010; Lu, Jensen, Ramani, & Helpern, 2006), Diffusion Spectrum Imaging (DSI; Wedeen, Hagmann, Tseng, Reese, & Weisskoff, 2005), q‐ball imaging (QBI; Tuch, Reese, Wiegell, & Wedeen, 2003), and high angular resolution diffusion imaging (HARDI; Tuch et al., 2002; Webster & Descoteaux, 2015) have been developed to index white matter alterations (Chanraud, Zahr, Sullivan, & Pfefferbaum, 2010; Koirala et al., 2019), neural density (Koirala, Perdue, Su, Grigorenko, & Landi, 2021; Zhang, Schneider, Wheeler‐Kingshott, & Alexander, 2012), and estimates of white matter fiber tracts (Jbabdi & Johansen‐Berg, 2011; Koirala et al., 2016). However, as diffusion and the generated indices (e.g., Fractional anisotropy, mean diffusivity) are influenced by several anatomical factors such as fiber arrangements, degree of myelination, cell membranes, microtubules, and axonal integrity, caution must be taken for interpreting the changes associated with these indices (Alba‐Ferrara & de Erausquin, 2013; Beaulieu, 2009; Jones & Cercignani, 2010; Jones, Knosche, & Turner, 2013).

Each of these techniques has expanded our understanding of structure–function relationships and various neurological and psychological disorders in both adults and children. Over the last decade, dMRI techniques have been increasingly employed in studies of neurodevelopmental disorders (e.g., attention deficit hyperactivity disorder (ADHD) (Ameis et al., 2016; van Ewijk, Heslenfeld, Zwiers, Buitelaar, & Oosterlaan, 2012; Wu et al., 2017), dyslexia (Vandermosten, Boets, Wouters, & Ghesquiere, 2012; Wang et al., 2017; Yeatman, Dougherty, Ben‐Shachar, & Wandell, 2012), autism spectrum disorder (ASD; Andrews et al., 2019; Ismail et al., 2016; Travers et al., 2012), and obsessive–compulsive disorder (OCD; Gruner et al., 2012; Jayarajan et al., 2012; Silk, Chen, Seal, & Vance, 2013, among others). While such studies have been helpful in elucidating differences in microstructural properties in these populations, there has been no systematic evaluation of how variable signal to noise ratio (SNR) and/or the presence of artifact relates to any of the dependent variables of interest, despite the fact that reports consistently note more motion and artifacts among children, especially those with neurodevelopmental disorders (Afacan et al., 2016; Dosenbach et al., 2017; Greene et al., 2018).

In general, compared to other structural imaging techniques, dMRI datasets are often characterized by low signal to noise ratio (SNR; Chilla, Tan, Xu, & Poh, 2015; Polders et al., 2011) and are frequently corrupted by multiple artifacts typically originating from eddy currents, insufficient fat‐suppression, B0 inhomogeneity, physiologically related factors (e.g., respiratory motion, cardiac pulsation, participant motion), and Gibbs ringing (Le Bihan, Poupon, Amadon, & Lethimonnier, 2006; Perrone et al., 2015; Pierpaoli, 2012; Tournier, Mori, & Leemans, 2011). Understandably, these artifacts, and particularly motion‐related artifacts, are even more significant for the studies involving young children and individuals with neurological developmental disorders (Tamnes, Roalf, Goddings, & Lebel, 2018; Theys, Wouters, & Ghesquiere, 2014). Moreover, it is increasingly common for studies to acquire data from multiple scanners and centers to increase sample size and sample diversity [i.e., Human Connectome Project (Van Essen et al., 2012), Adolescent Brain Cognitive Development study (Casey et al., 2018), etc.]. Pooling data collected from different centers brings new challenges of harmonizing data collection techniques, and data quality from different scanners may be differentially impacted by variations in hardware, scanning sequences, environmental factors, and human factors, as well as other sources of variability (Fortin et al., 2017; Mirzaalian et al., 2016; Ning et al., 2020; Pinto et al., 2020). This is especially relevant when pooling data from existing datasets in which data collection was not optimized for combination with data from other sources. Thus, accounting for cross‐scanner differences in image quality when analyzing data collected at multiple centers is a critical step in data analysis.

Over the years, there have been several advances in the attempt to overcome scanning artifacts with different protocols (Bammer, Holdsworth, Veldhuis, & Skare, 2009; Nana, Zhao, & Hu, 2008; Nolte, Finsterbusch, & Frahm, 2000), motion correction techniques (Alhamud, Taylor, Laughton, van der Kouwe, & Meintjes, 2015; Chan et al., 2014; Truong, Chen, & Song, 2011), and better processing algorithms (Andersson et al., 2017; Bastiani et al., 2017; Li et al., 2014). However, there is no consensus on how to compensate for differences in acquisition protocol, artifact handling, data quality control, reconstruction algorithm, visualization approach, or quantitative analysis methodology. Moreover, previous studies have shown that diffusion measures such as the apparent diffusion coefficient (ADC) and fractional anisotropy (FA) vary substantially with different vendors, coil systems, imagers, field strengths, and sequence parameters (Helmer et al., 2016; Sasaki et al., 2008; Schmeel, 2019). The inter‐subject, inter‐session, inter‐site, and between‐visits variability were also shown to be significant for these measures with varying degree of significance (Huo et al., 2016; Veenith et al., 2013). In addition, the assumptions made during the acquisition (e.g., perfect field homogeneity, infinitely fast gradient changes, perfectly shaped RF pulses) and processing (e.g., complete correction of susceptibility‐induced distortion and eddy current effects, perfect alignment of regions of interests during registration) of DWI data have many pitfalls, and these shortcomings have been shown to effect the computed diffusion measures (Baliyan, Das, Sharma, & Gupta, 2016; Hrabe, Kaur, & Guilfoyle, 2007; Jones & Cercignani, 2010). Even though it is generally accepted that image quality has an impact on the computed measures (Bastin, Armitage, & Marshall, 1998; Soares, Marques, Alves, & Sousa, 2013) and on diffusion measures specifically (Baum et al., 2018; Landman et al., 2007; Ling et al., 2012), research investigating the full range of diffusion measures and the extent of the impact of image quality is still sparse. Here, we investigate (a) several data quality metrics (the ones most relevant to diffusion imaging data: different components of motion during scanning—absolute, translational, rotational, and relative—as well as contrast to noise ratio and signal drop in different slices, or the number of outlier slices) obtained from children of age range of 5–17 years with varying disability status, and (b) the association of these quality metrics with computed diffusion (fractional anisotropy and mean diffusivity) and tractography (tract density and length) measures obtained from a multicenter dataset with different diffusion modalities and magnetic strength.

METHODS

Data acquisition

To assess the relationships between data quality and diffusion measures, we analyzed 708 dMRI data sets of 691 participants who were between 5 and 17 years of age (17 participants in the Italy dataset were scanned twice with two different protocols). Data were obtained from six different centers, using four different models of scanner, five different scanning sequences, two different magnetic strengths (1.5 Tesla [T] and 3 T), and containing two different diffusion models (DTI and DKI). dMRI datasets were obtained from Haskins Laboratories (HL), Boston Children Hospital (BCH), IRCCS Eugenio Medea Italy (ITA), and the Healthy Brain Network Biobank (L. M. Alexander et al., 2017), which includes data from three centers: CitiGroup Cornell Brain Imaging Center (CBIC), Rutgers University Brain Imaging Center (RU) and Staten Island (SI). Sample sizes, demographic information, and data acquisition parameters are presented in Table 1 and Figure 1. All participants included in the study were participating in a multisite study as part of the Florida Learning Disabilities Research Center project VI: Imaging genetics in SRD: Mega‐ and meta‐analyses. The data acquisition procedure was approved by the Chesapeake Institutional Review Board (https://www.chesapeakeirb.com/) for the Healthy Brain Network dataset and respective ethics commission for each center. Prior to acquiring the data, written assent obtained from the participant, and written informed consent was obtained from their legal guardians (L. M. Alexander et al., 2017).

TABLE 1

Details of scanning parameters from all centers

Centers	Number of participants	Mean age (years)	Scanner (field strength)	Diffusion model	Diffusion directions	b values (s/mm2)	Voxel size (mm)
CBIC	228	9.98 ± 2.96	Siemens Prisma (3 T)	DKI	64	0, 1,000, 2,000	2.0 × 2.0 × 2.0
RU	233	9.36 ± 2.86	Siemens TrioTim (3 T)	DKI	64	0, 1,000, 2,000	2.0 × 2.0 × 2.0
SI	117	10.96 ± 3.10	Siemens Avanto (1.5 T)	DKI	64	0, 1,000, 2,000	2.0 × 2.0 × 2.0
HL	37	6.91 ± 0.89	Siemens TrioTim (3 T)	DTI	32	0, 1,000	2.0 × 2.0 × 3.0
BCH	32	5.39 ± 0.27	Siemens TrioTim (3 T)	DTI	30	0, 1,000	2.0 × 2.0 × 2.0
ITA_1	34	13.19 ± 1.94	Philips Achieva (3 T)	DTI	32	0, 300, 1,100	1.7 × 1.7 × 2.0
ITA_2	27	12.90 ± 1.79	Philips Achieva (3 T)	DTI	32	0, 1,100, 2,500	1.7 × 1.7 × 2.0

Note: Here, ITA_1 and ITA_2 is separately shown to highlight the data acquisition difference in scanning parameters (b values) used for subanalysis.

FIGURE 1

Age and sex distribution of all the participants. All participants were selected without limitations on sex or gender, race, or ethnicity, or age other than as scientifically justified. Here M and F indicates male and female and the black dash shows the median line

Reproducibility and data availability

All the raw data (Neuroimaging) used in the study from sites CBIC, RU, and SI are freely available via Healthy Brain Network Biobank (https://childmind.org/center/healthy-brain-network/). The phenotypical data used could be available upon request. The toolbox used in the study (FSL) is an open‐access toolbox that can be downloaded from https://fsl.fmrib.ox.ac.uk/fsldownloads_registration. For specific processing protocols used in the study, please refer to Section 2 and contact the corresponding author for any further queries. Scripts used to conduct the analyses reported in this article have been made available on the Open Science Framework (OSF) at https://osf.io/9s27f.

Data analysis

The MRI scans primarily obtained in DICOM format were first converted into Brain Imaging Data Structure (BIDS) format using the dcm2bids toolbox available freely at https://github.com/cbedetti/Dcm2Bids (Gorgolewski et al., 2016). Data quality metrics were computed for all participants using an automated quality control framework (named QUAD—QUality Assessment for DMRI; Bastiani et al., 2019) available as an open‐source toolbox FSL (ver. 6.0.3). The diffusion data quality metrics (dMRIqc) include average‐absolute motion (AAM), average‐relative motion (ARM), average‐translational motion (ATM), and average‐rotational motion (AOM) motion, contrast to noise ratio (CNR), and the number of outlier slices (outlier_dwi). Average absolute motion (with respect to a reference volume) and relative motion (with respect to the previous volume) were calculated as the average voxel displacement across all voxels within a brain mask summarizing both translations and rotations at each voxel (Bastiani et al., 2019). The number of outlier slices indexed how many single slices got distorted during data acquisition because of participant motion that caused signal drop (Andersson & Sotiropoulos, 2016). Here, we computed the total percentage of slices that were distorted within a given volume for each participant.

To obtain the diffusion and tract measures, the images were preprocessed using inbuilt functionality in the open‐source toolbox FSL (ver. 6.0.3) described in detail elsewhere (Jenkinson, Bannister, Brady, & Smith, 2002; Jenkinson & Smith, 2001). In brief, data collected with reversed phase‐encode blips (except for site BCH) were used for estimating and correcting the susceptibility‐induced artifacts using topup (Andersson, Skare, & Ashburner, 2003; Smith et al., 2004). Corrections of motion artifacts (eddy currents and head movements) were performed using the eddy toolbox (Andersson & Sotiropoulos, 2016), and individual masks were generated for each brain using the Brain Extraction Toolkit (BET; Smith, 2002) to isolate the brain from the skull and diffusion tensor modeling for obtaining diffusion measures. Fractional anisotropy (FA) and Mean diffusivity (MD) were calculated using the FDT toolbox. In addition, the distribution of crossing fibers was estimated using BEDPOSTX (T. Behrens, Berg, Jbabdi, Rushworth, & Woolrich, 2007; T. E. Behrens et al., 2003), and the probability of major (f1) and secondary (f2) fiber directions was calculated (Koirala et al., 2017; Koirala et al., 2019). The obtained crossing fiber modeled diffusion data were further processed using the automatic tractography scheme using XTRACT toolbox in FSL (Warrington et al., 2020). Here, tractography masks are defined in standard space; these masks were then warped to each participant's native space using the participant‐specific, nonlinear warp fields; and probabilistic tractography was performed in the participant's native space. The resultant tract was then stored in standard space and overlaid on the FSL_HCP1065 FA atlas (Warrington et al., 2020). Finally, tract density, tract length, FA, and MD values were extracted from each of the 36 major tracts as detailed in Table 2.

TABLE 2

All tracts where the diffusion parameters—tract—density, length, fractional anisotropy (FA), and mean diffusivity (MD) values were extracted

Tract	Abbreviation
Anterior commissure	AC
Arcuate fasciculus (left and right)	Left AF, right AF
Acoustic radiation	AR
Anterior thalamic radiation (left and right)	Left ATR, right ATR
Dorsal cingulum (left and right)	Left CBD, right CBD
Peri‐genual cingulum (left and right)	Left CBP, right CBP
Temporal cingulum (left and right)	Left CBT, right CBT
Corticospinal tract (left and right)	Left CST, right CST
Forceps major	FMA
Forceps minor	FMI
Fornix (left and right)	Left FX, right FX
Inferior longitudinal fasciculus (left and right)	Left ILF, right ILF
Inferior fronto‐occipital fasciculus (left and right)	Left IFO, right IFO
Middle longitudinal fasciculus (left and right)	Left MdLF, right MdLF
Optic radiation (left and right)	Left OR, right OR
Superior longitudinal fasciculus 1 (left and right)	Left SLF1, right SLF1
Superior longitudinal fasciculus 2 (left and right)	Left SLF2, right SLF2
Superior longitudinal fasciculus 3 (left and right)	Left SLF3, right SLF3
Uncinate fasciculus (left and right)	Left UF, right UF
Vertical occipital fasciculus (left and right)	Left VOF, right VOF

Using the obtained parameters, statistical analyses were conducted to detect differences between these data quality metrics across centers and brain regions (tracts) and to determine the relationship between these data quality metrics and diffusion measures, as detailed in Section 2.4. In recent years, different harmonization algorithms have been proposed for removing unwanted inter‐site variability from multi‐site diffusion data (Fortin et al., 2017; Mirzaalian et al., 2016; Ning et al., 2020; Pinto et al., 2020). These algorithms reduce/remove the systematic differences between scanner manufacturers, field strength, and other scanner characteristics that systematically affect the diffusion images and introduce inter‐scanner variation. A recent study (Fortin et al., 2017) attempted to address this issue by comparing five statistical harmonization techniques for diffusion data: global scaling, functional normalization (Fortin et al., 2014), removal of artificial voxel effect by linear regression (RAVEL; Fortin et al., 2016), surrogate variable analysis (SVA; Leek & Storey, 2007), and ComBat (Johnson, Li, & Rabinovic, 2007). Of these, the ComBat model, which uses an empirical Bayes (EB) framework to improve the variance of the parameter estimates, was shown to be the most effective in removing unwanted variation induced by site in diffusion data. Accordingly, to evaluate whether applying a harmonization technique in our multi‐site dataset would affect the results of our study, we harmonized our data using the ComBat model. To do this, we applied the Matlab‐based algorithm (open access availability at https://github.com/Jfortin1/ComBatHarmonization) to harmonize the obtained diffusion measures—tract density, tract length, FA, and MD—across sites. A second set of statistical analyses was run using the data harmonized via ComBat.

Statistical analysis

Four sets of analyses were conducted in R (v. 3.6.1; R Core Team, 2019) to examine the relationships between six data quality metrics (AAM, ARM, ATM, AOM, CNR, and outlier_dwi) and four diffusion measures (tract density, tract length, FA, and MD) within and across centers. To examine the effects of interest for each set of analyses, post hoc contrasts were applied to the parameters estimated by the statistical models.

Analysis 1

To determine whether data quality metrics differed by center, six multiple regressions were conducted—one for each data quality metric—to compute effects of center on that metric (the dependent variable). To explore differences in these metrics between centers in more detail, pairwise comparisons were conducted between centers for each metric, with the Tukey method used to control family‐wise error rate when the main effect of center was significant.

Analysis 2

To determine whether diffusion measures differed by center and across brain regions, four linear mixed‐effects regression models (Baayen, Davidson, & Bates, 2008) were fit—one for each diffusion measure—to compute effects of center, tract, and their interaction on that diffusion measure (the dependent variable). Contrasts measured the extent to which centers differed in diffusion measures at each of the fiber tracts obtained.

Analysis 3

To determine the relationship between data quality metrics and diffusion measures as well as any modulation of this relationship by center and brain region, 24 mixed‐effects models were fit—one for each combination of data quality metric (6) and diffusion measure (4)—to compute effects of center, tract, data quality metric, and their interactions on that diffusion measure (the dependent variable). Contrasts examined the relationship between each pair of data quality and diffusion measures at each of the fiber tracts obtained.

Analysis 4

To determine the impact of different scanner models, magnetic strength, diffusion model, and sequence on the relationship between data quality metrics and diffusion measures, we performed contrasts on the statistical models fitted for Analysis 3 (no new models were fit for this analysis). Each contrast compared two centers (or two groups of centers), examining the extent to which the (groups of) centers differed on the relationship between data quality metrics and diffusion measures—both across all tracts, separately at each individual tract. These contrasts examined the effects of several variables:We note that although the centers in each contrast were selected to maximize the relevant difference while holding other variables as constant as possible, all between‐center contrasts (A, B, and C) likely captured the joint influences of several variables on diffusion metrics.

Scanner models (Siemens Prisma vs. TrioTim): By comparing CBIC (N = 228, Siemens Prisma) with RU (N = 233, Siemens TrioTim).

Magnetic strength (3 vs. 1.5 T): By comparing CBIC (N = 228, 3 T) and RU (N = 233, 3 T) with SI (N = 117, 1.5 T).

dMRI models (DKI vs. DTI): By comparing CBIC (N = 228, DKI) with HL (N = 37, DTI) and BCH (N = 32, DTI).

dMRI sequence (higher vs. lower b values ): By comparing ITA_1 and ITA_2, where ITA_2 was a partially overlapping subset of participants (overlapping n = 17) who were scanned during the same acquisition period but using higher b values acquisition sequence (Table 1).

All analyses included fixed effects of age and sex, both as main effects and in interactions with all other fixed effects. Including these effects in every model ensured that wherever the dependent variable covaried with age and/or sex, that variance was not attributed to critical predictors (e.g., effects of center or tract or data quality metric). Also, for all analyses, contrast weightings for the effect of center were selected to weight each center in proportion to the number of participants tested at that center (a form of weighted effects coding; cf. [te Grotenhuis et al., 2017]). This means that, for example, when evaluating the main effect of tract, the data from each participant was weighted the same regardless of whether they were tested at center RU (n = 233) or center HL (n = 32)—which in turn means that, overall, the data from RU was weighted 7.3 times as much as from HL (because 7.3 times as many participants were tested there).

For Analysis 1, which used multiple regressions, significance tests for each factor were conducted via model comparison using the ANOVA function. For Analysis 2 and Analysis 3, which used mixed‐effects models, significance tests were conducted using the contestMD function in the lmerTest package (v. 3.1‐2; Kuznetsova, Brockhoff, & Christensen, 2017), using the Satterthwaite method to approximate denominator degrees of freedom. (These methods of significance testing were used because they preserved the contrast weights for each center during significance testing.) The mixed‐effects models themselves were conducted using the lme4 package (v. 1.1–21; Bates, Machler, Bolker, & Walker, 2015), employing a two‐step model fitting strategy. First, a model was fit with a maximal random effects structure by adding a random intercept for participants, as well as all within‐participant random slopes and their interactions. (Correlations between random slopes were not added because this would have caused the number of random effects parameters to equal the number of observations, making the model unidentifiable. Also, as 674/691 participants (97.5%) were tested at one center each—that is, the entire sample except for the 17 participants tested at the Italy site with two b values acquisition sequences—all scans were treated as independent for analytic purposes.) If this model did not converge, all random slopes accounting for less than 1% of the random variance were removed simultaneously (Bates et al., 2015), which always resulted in convergence. Contrasts were computed using the emmeans package (v. 1.4; Lenth & Love, 2018). To control Type I error rate, corrections for multiple comparisons were applied. Within each analysis, except as noted below, the Benjamini–Yekutieli method (Benjamini & Yekutieli, 2001) was used to set the false discovery rate (FDR) at 0.05 across all diffusion measures, data quality metrics, and tracts. This correction was applied separately to each effect of interest because each effect represents a different hypothesis. For example, we conducted 24 tests of the hypothesis that there is a main effect of data quality metric on diffusion measures (one test for each combination of six metrics and four measures), so we applied a correction for 24 tests to those p values. Separately, we conducted 24 tests of the hypothesis that the effect of data quality metric on diffusion measures interacts with center, so we applied a correction for 24 tests to the p values of that two‐way interaction. (For each hypothesis test, the number of comparisons corrected for that test is included in the Supporting Information.)

Testing the effect of data harmonization

We then ran a second set of analyses on the ComBat‐harmonized metrics, including the analyses examining relationships between data quality metrics and diffusion metrics at different centers and tracts. These analyses were identical to our normal analyses—including effects of age and sex on top of the harmonization—except that we also removed the data from center ITA2 to avoid harmonizing data from a subset of participants twice.

RESULTS

Analysis 1

We found that four data quality metrics (all except for AAM and ATM) significantly differed between centers, as indicated by a significant main effect of center in four of the six models. Test statistics for each data quality metric are given in Table 3. To explore differences in these metrics between centers in more detail, pairwise comparisons were conducted between centers for each metric, with the Tukey method used to control family‐wise error rate when the main effect of center was significant. These pairwise comparisons are shown in Figure 2.

TABLE 3

Details of the results obtained from Analysis 1—To determine whether data quality metrics differed by center

Image quality metric	Between centers
Average relative motion (ARM)	F (6,687)  = 30.2, p < .001
Contrast‐to‐noise ratio (CNR)	F (6,687)  = 25.2, p < .001
Outlier slices (outlier_dwi)	F (6,687)  = 13.8, p < .001
Average rotational motion (AOM)	F (6,687)  = 4.7, p < .001
Average absolute motion (AAM)	F (6,687) = 0.8, p = 1.00
Average translation motion (ATM)	F (6,687) = 0.5, p = 1.00

Note: The Benjamini–Yekutieli method was used to control FDR across image quality metrics; adjusted p values are reported. Significant results are shown in bold.

FIGURE 2

Data quality metrics for all centers. For metrics that significantly differed across centers, letters above each center show the results of Tukey‐corrected pairwise comparisons: Centers that do not share a letter significantly differed on the metric after regressing out age and sex. For example, center CBIC significantly differed from centers SI and ITA2, but not any other center, on contrast‐to‐noise ratio

Analysis 2

We found that all analyzed diffusion measures significantly differed between centers and between tracts, and that these effects interacted. Test statistics for each diffusion measure are presented in Table 4. Descriptively, centers differed most strongly on MD, followed by tract density, tract length, and FA (Figure 3), whereas tracts differed most strongly on tract length and FA, followed by MD and tract density. All 36 tracts significantly differed across centers on FA, MD, and tract length; and 11/36 tracts significantly differed across centers on tract density (Figure 4).

TABLE 4

Details of the results obtained from Analysis 2—To determine whether diffusion measures differed by center and across brain regions

Diffusion measure	Between centers	Between tracts	Center × tract
FA	F (6,687)  = 17.7, p < .001	F (35,1,352)  = 136.4, p < .001	F (210,1,352)  = 35.7, p < .001
MD	F (6,687)  = 595.0, p < .001	F (35,1,604)  = 71.6, p < .001	F (210,1,604)  = 48.0, p < .001
Tract density	F (6,687)  = 188.3, p < .001	F (35,2,746)  = 53.3, p < .001	F (210,2,746)  = 99.3, p < .001
Tract length	F (6,687)  = 64.0, p < .001	F (35,1,159)  = 170.6, p < .001	F (210,1,159)  = 37.9, p < .001

Note: The Benjamini–Yekutieli method was used to control FDR across diffusion measures; adjusted p values are reported. Significant results are shown in bold.

FIGURE 3

Diffusion measures for all centers. For measures that significantly differed across centers, letters above each center show the results of Tukey‐corrected pairwise comparisons: Centers that do not share a letter significantly differed on the measure after regressing out age and sex. For example, center CBIC significantly differed from center SI, but not any other center, on tract length. As analyses controlled for participant age and sex, centers that were outliers in their distribution of participant ages (e.g., center BCH) may show fewer significant differences with other centers here than their means would suggest

FIGURE 4

All tracts that differ significantly across centers on Mean diffusivity (MD) (36/36), fractional anisotropy (36/36), tract length (36/36), and tract density (11/36). Each line indicates a tract with the height of it showing strength. Please refer to Supporting Information for the details of all tracts shown here with the statistical values

Analysis 3

Associations between diffusion measures and data quality metrics varied across tracts to different degrees at different centers. This was indicated by significant three‐way interactions between center, tract, and data quality metrics for 13/24 (54%) combinations of diffusion measures and data quality metrics (Figure 5, bottom row). Relationships with MD were the most influenced, showing significant three‐way interactions with center and tract on 5/6 data quality metrics, followed by tract density (4/6 data quality metrics). The least impacted relationships were with FA and tract length, for which center and tract jointly modulated relationships with 2/6 data quality metrics each (both average CNR and ARM). Relatedly, we found that relationships with ARM were significantly modulated by center and tract for 4/4 diffusion measures, and relationships with CNR were significantly modulated by center and tract for 3/4 diffusion measures (all except tract density).

FIGURE 5

Figure visualizing Analysis 3 effects of data quality metrics on diffusion measures. Each row corresponds to a different type of effect (main effect of a data quality metric in the top row, two‐way interactions with center or tract in the middle rows, three‐way interaction with both center and tract in the bottom row). Each column corresponds to a different diffusion measure. Within each panel, each colored bar corresponds to a different data quality metric. The Y‐axis represents the F‐value of the given effect with the given data quality metric on the given diffusion measure. Effects that were statistically significant after correction for multiple comparisons are shaded with diagonal slashes. (Bar heights are directly comparable within a row but are not comparable between rows due to differences in degrees of freedom for different effects.) As an example of how to interpret the figure, the prevalence of many shaded bars in the MD column means that multiple significant relationships were observed between data quality metrics and MD, whereas the prevalence of many yellow shaded bars throughout the figure means that CNR was significantly related to multiple diffusion measures in different ways. Please refer to Supporting Information for the details of all tracts shown here with the statistical values. AAM, average absolute motion; AOM, average rotational motion; ARM, average relative motion; ATM, average translational motion; CNR, contrast to noise ratio; FA, fractional anisotropy; MD, mean diffusivity

Although the presence of three‐way interactions qualifies the interpretation of two‐way interactions, we note for completeness that for MD and tract length, relationships with several data quality metrics significantly differed across both centers and tracts. For the remaining diffusion measures, however, significant relationships with one or more data quality metrics were observed either across centers (FA) or across tracts (tract density; see Figure 5, middle rows).

A small number of significant relationships (11/864, after correcting for multiple comparisons) were observed between diffusion measures and image quality metrics within individual tracts. Positive relationships were observed between CNR and MD in right optic radiation, left superior longitudinal fasciculus 1, and bilateral uncinate fasciculus; and between CNR and tract length in left corticospinal tract, left inferior fronto‐occipital fasciculus, right superior longitudinal fasciculus 1, left superior longitudinal fasciculus 2, and bilateral uncinate fasciculus; in addition, a negative relationship was observed between the number of outlier slices and MD in the forceps major. The results obtained from these analyses provide further evidence that the relationships between data quality metrics and diffusion measures are irregularly modulated across different tracts (even if they are only significant at a few individual tracts) and vary significantly across centers.

Analysis 4

Scanner models (Siemens Prisma vs. TrioTim)

When comparing centers differing in the scanner manufacturer (CBIC vs. RU), we found several significant differences in the relationships between diffusion measures and data quality metrics, both across tracts and at individual tracts. Modulation across tracts was observed for the relationship between MD and image quality metrics (ARM, ATT, and CNR), as well as the relationships between FA and CNR, and between tract length and ARM (Table 5). Moreover, the relationships between the metrics varied between tracts, with the largest number of significant effects (22) associated with MD, then tract length (6) and FA (5) [and none with tract density]. All tracts at which a significant difference was observed between centers in the relationship between diffusion measures and data quality metrics are shown in Figure 6.

TABLE 5

Details of significant (p < .05) results from Analysis 4—To determine the impact of different scanner models, magnetic strength, diffusion model, and sequence to the relationship of data quality metrics and diffusion measures

Scanner models
Diffusion measure	Data quality metric	B ± SE
FA	CNR	−0.02 ± 0.004
MD	ARM	1E−5 ± 3E−6
MD	ATM	1E−5 ± 3E−6
MD	CNR	1E−5 ± 5E−6
Tract length	ARM	1.64 ± 0.513

Note: The Benjamini–Yekutieli method was used to control FDR across diffusion measures. Here, B indicates the Beta coefficient representing the mean effect of the predictor (data quality metrics) on the dependent variable (DV; diffusion measures) for one center minus the effect of the predictor on the DV for another center. The overall contrast for this analysis was only significant for scanner models. See full results provided at https://osf.io/9s27f.

Abbreviations: ARM, average relative motion; ATM, average translation motion; CNR, contrast to noise ratio; FA, fractional anisotropy; MD, mean diffusivity; SE, standard error.

FIGURE 6

All tracts that showed statistically significant differences between the centers in the effects of each data quality metric in each tract on each diffusion measures. Please refer to Supporting Information for the details of all tracts shown here with the statistical values. Here, B indicates the Beta coefficient representing the mean effect of the predictor (data quality metrics) on the dependent variable (DV) (diffusion measures) for one (group of) center(s) minus the effect of the predictor on the DV for another (group of) center(s). AAM, average absolute motion; AOM, average rotational motion; ARM, average relative motion; ATM, average translational motion; CNR, Contrast to noise ratio; FA, fractional anisotropy; MD, mean diffusivity

Magnetic strength (3 vs. 1.5 T)

When comparing centers differing in magnet strength (CBIC and RU vs. SI), no significant differences were observed when averaging across tracts. Comparing centers at each tract, all three significant differences were in the relationship between ARM and various diffusion metrics in the temporal cingulum, with stronger relationships observed for the 1.5 T scanner. The details of all tracts showing significant differences are shown in Figure 6. (Note that the centers in this contrast also varied on scanner model, as shown in Table 1.)

dMRI models (DKI vs. DTI)

When comparing centers that differed in the diffusion models (DKI and DTI) used to acquire data (CBIC vs. HL and BCH), no significant differences were observed when averaging across tracts. Comparing centers at each tract, one significant difference was observed in the relationship between outlier_dwi and MD at FMA, with a stronger relationship observed for the DTI model. (Note that the centers in this contrast also varied on scanner model, the number of diffusion directions, b values, and voxel size, as shown in Table 1.)

dMRI sequence (higher vs. lower b values)

For this dataset acquired during the same session with two different b values (ITA1 vs. ITA2), we did not observe any significant differences when averaging across tracts, but significant differences did emerge in the relationships between data quality metrics and two diffusion measures—MD and tract density—at a small number of tracts (Figure 6). These effects were not always more prominent for the higher b values sequence (probably as one might assume) but were rather dependent on the data quality metric and diffusion measures.

The effect of harmonization

The analyses with ComBat‐harmonized data revealed that the harmonization did not eliminate or even reduce the number of significant relationships observed between the six data quality metrics and the four diffusion measures. In Analysis 3, those relationships

were significant (main effect) for 4/24 combinations of metrics (vs. 2/24 for non‐harmonized data);

significantly varied by center (two‐way interaction) for 0/24 combinations of metrics (vs. 5/24 for non‐harmonized data);

significantly varied by tract (two‐way interaction) for 5/24 combinations of metrics (vs. 6/24 for non‐harmonized data); and

significantly varied jointly by center and tract (three‐way interaction) for 17/24 combinations of metrics (vs. 13/24 for non‐harmonized data)

Regarding tract‐by‐tract analyses, of all combinations of 36 tracts, six data quality metrics, and four diffusion metrics, the relationship between metrics was significant at an individual tract for 12/864 combinations (vs. 11/864 for non‐harmonized data).

In general, ComBat harmonization did appear to reduce between‐center differences in diffusion metrics (the main effect of center and the interaction between center and tract were each significant for 0/4 metrics in Analysis 2, vs. 4/4 for non‐harmonized data). However, perhaps surprisingly, the number of significant interactions with age and sex in all analyses was roughly comparable for harmonized and non‐harmonized data.

All the statistical results from the Main as well as from ComBat analysis have been provided as Supporting Information in the OSF platform (https://osf.io/9s27f).

DISCUSSION

In this study, involving a large multisite dataset of typically developing children and children with neurodevelopmental disorders, we demonstrated that the diffusion imaging data quality, assessed via multiple metrics, differed significantly between data collection sites. Furthermore, the associations between these data quality indices and the obtained diffusion measures differed non‐uniformly across sites and brain areas. Analyzing different components of data acquisition, we found significant modulation of the relationship between data quality metric and diffusion measures were significantly influenced by scanner models, magnetic strength, diffusion models, and data collection sequences.

Over the last two decades, dMRI has been established as an effective technique for understanding various neurodevelopmental conditions and disorders and psychological states and processes. Several indices have been proposed encompassing a wide range of microstructural properties including anisotropy, diffusivity, fiber orientation, and neurite features. These indices are known to be influenced by several white matter properties such as fiber arrangements, degree of myelination, and axonal integrity (Alba‐Ferrara & de Erausquin, 2013; Jones et al., 2013). Moreover, the contributions of these white matter properties to computed diffusion measures are non‐linear and not well known. To give an example, in a recent study using animal models treated with different drugs to trigger demyelination and axonal damage, the FA values obtained could not distinguish between a demyelination group and those with both demyelination and axonal damage (Boretius et al., 2012). Similar concerns have been raised for crossing fibers (around 90% of white matter voxels contain crossing fibers) and q‐space (diffusion space), which have been shown to have a remarkable impact on fiber tract estimation and anisotropy analysis (Jeurissen, Leemans, Tournier, Jones, & Sijbers, 2013; Wilkins, Lee, Gajawelli, Law, & Lepore, 2015). The susceptibility of dMRI data to various artifacts (e.g., eddy currents, echo‐planar distortions, rotation errors, partial volume effects, and scanner artifacts) has been studied to better understand how these artifacts impact the estimation of eigenvalues and eigenvectors, which in turn influence the accuracy of computed anisotropic measures and fiber tracking schemes (Anderson, 2001; Skare, Li, Nordell, & Ingvar, 2000). Furthermore, one of the largest confounds for studies in children and adolescents—head motion—has been shown to have a significant impact on results obtained in group comparison studies (Yendiki, Koldewyn, Kakunoori, Kanwisher, & Fischl, 2014) or in studies investigating the relationship to diffusion measures (Baum et al., 2018). Even though several studies have highlighted these aspects and proposed new methods (Li et al., 2014; Oguz et al., 2014), the consideration of data quality for the analysis and a transparent, standardized estimate of quality assurance is still rare. Most of the studies suggest either removing the artifacts during preprocessing steps or manually performing quality control as a good practice for controlling for this issue. In this study, we provide evidence that even after these standard preprocessing steps for artifacts correction, the relationships between these data quality measures and computed diffusion metrics are still widespread across the brain and vary non‐uniformly across tracts and data collection centers. This result has potentially important implications for the inferences made by the clinical and neurobiological studies that use these metrics.

The results from Analysis 1 and Analysis 2 illustrated that data quality metrics (except for absolute and translational motion) and computed diffusion measures varied across centers. This provides one possible scientific explanation for differences in findings and replication problems for studies conducted across different centers. Moreover, we observed that the impact of these data quality metrics on diffusion measures varies across centers and across different tracts in the brain (Analysis 3). The tracts which showed the strongest variation across centers for different diffusion measures are commonly studied in clinical populations and in relation to various neuropsychological parameters. In addition, we observed that MD, tract density, tract length, and FA, which are some of the widely used parameters in neurodevelopmental studies, were the parameters highly influenced by the data quality metrics. However, to the best of our knowledge, no prior neurodevelopmental study has accounted for these data quality metrics in the statistical models. A failure to do so might lead to unaccounted impacts on suboptimal tensor estimation, erroneous partial volume effect computation, or estimation of boundary threshold. Accordingly, we recommend computing these quality metrics and including them in analyses when reporting delicate microstructural changes.

In Analysis 4, the four subanalyses completed to quantify these differences across different scanner models, magnetic strength, diffusion models, and q‐space sampling or b values, revealed significant modulation of the relation between quality metrics and diffusion measures by these factors. For these comparisons, MD was the measure that showed the largest number of significant relationships with data quality measures, and ARM was the data quality metric that showed the largest number of significant relationships with diffusion measures. More importantly, the impact of data quality varied (e.g., direction of association,: positive or negative and strength of the association) across different tracts and different diffusion measures, making it impossible to process out one or more effects. For example, the effect of CNR in FA values is higher in most tracts for Siemens TrimTrio; however, the effect of the same in MD is higher for Siemens Prisma. In addition, the strength of the relationship between the FA, MD, tract length and tract density is not uniform for each of the data quality metrics and for each tract (refer to Figure 6 for visualization of these examples). This variation further highlights the importance of using the quality metrics in statistical analysis while computing these diffusion measures to capture this unseen influence. Moreover, for multisite data, this might be a promising way to gauge differences in the datasets and evaluate them for harmonization.

After harmonizing our data across centers using the ComBat algorithm, we continued to observe significant relationships between data quality metrics and diffusion metrics, both directly and varying across centers and tracts. Of these effects, the total number of significant effects was the same with vs. without harmonization (26/96 significant). This further supports the necessity of considering data quality metrics in data analysis even after harmonization.

In recent years, there have been several studies which have highlighted the impact of quality assurance of dMRI data on computed measures (Maximov, Alnaes, & Westlye, 2019; Roalf et al., 2016) and have suggested various optimization methods (Liu et al., 2019) and alternative diffusion measures less sensitive to the data quality (Ozcan, 2010). With this study, we presented quantitative sets of data quality metrics (that can be straightforwardly obtained using standard toolboxes) and comprehensive details of their impact on the computed diffusion measures. We suggest that moving forward, these quality metrics should be considered in data analysis, particularly for multisite studies and those involving populations prone to movement. The dataset used for this analysis was aggregated post‐hoc; thus, data collected at each center consisted of a different set of participants. Moreover, participant diagnoses for different neurodevelopmental and neuropsychological disorders were not modeled during analysis as this was outside the scope of the study. Although it is possible that including diagnosis could change which specific relationships with data quality metrics were statistically significant, we think that the overall pattern of results—that is, the existence of significant relationships between these metrics and diffusion measures, often modulated by center and/or tract—would remain the same. Despite the aggregated nature of the dataset, we showed strong effects of center and influence of different components involved in data acquisition including scanner model, magnetic strength, and diffusion models (though some of these comparisons may reflect other between‐center variables as well). We conclude that this assessment of image quality in such a dataset is important as it most closely represents the types of large‐scale multisite studies that researchers currently draw upon for conducting large‐scale neuroimaging research. However, ongoing and future multisite studies can and should further optimize harmonization across sites through the use of concurrent acquisition models and sequence where possible, evaluate the difference with the use of phantoms beforehand, and use some of the harmonization techniques (Fortin et al., 2017; Mirzaalian et al., 2016; Ning et al., 2020) for closely related sequences.

CONCLUSION

We provided evidence that several dMRI data quality metrics differ between imaging data collection sites. More importantly, the relationship between these data quality metrics and diffusion measures are widespread across brain regions, though this too varied across center. The obtained brain areas and diffusion measures have been widely employed to discriminate neurotypical populations from individuals with neurological or developmental disorders. Hence, these data quality metrics should be considered for inclusion in dMRI analysis, which could be quantified using available toolboxes. Taking these metrics into account would help make data more comparable across testing sites, would provide more reliable findings, and might be a key in producing replicable studies for group comparison, identifying correlates of individual differences, and carrying out large‐scale multisite studies to identify neural markers for different neurodevelopmental conditions and cognitive functions.

CONFLICT OF INTEREST

All authors have no conflict of interest to declare.