Multisite Harmonization of Structural DTI Networks in Children: An A-CAP Study

The analysis of large, multisite neuroimaging datasets provides a promising means for robust characterization of brain networks that can reduce false positives and improve reproducibility. However, the use of different MRI scanners introduces variability to the data. Managing those sources of variability is increasingly important for the generation of accurate group-level inferences. ComBat is one of the most promising tools for multisite (multiscanner) harmonization of structural neuroimaging data, but no study has examined its application to graph theory metrics derived from the structural brain connectome. The present work evaluates the use of ComBat for multisite harmonization in the context of structural network analysis of diffusion-weighted scans from the Advancing Concussion Assessment in Pediatrics (A-CAP) study. Scans were acquired on six different scanners from 484 children aged 8.00–16.99 years [Mean = 12.37 ± 2.34 years; 289 (59.7%) Male] ~10 days following mild traumatic brain injury (n = 313) or orthopedic injury (n = 171). Whole brain deterministic diffusion tensor tractography was conducted and used to construct a 90 x 90 weighted (average fractional anisotropy) adjacency matrix for each scan. ComBat harmonization was applied separately at one of two different stages during data processing, either on the (i) weighted adjacency matrices (matrix harmonization) or (ii) global network metrics derived using unharmonized weighted adjacency matrices (parameter harmonization). Global network metrics based on unharmonized adjacency matrices and each harmonization approach were derived. Robust scanner effects were found for unharmonized metrics. Some scanner effects remained significant for matrix harmonized metrics, but effect sizes were less robust. Parameter harmonized metrics did not differ by scanner. Intraclass correlations (ICC) indicated good to excellent within-scanner consistency between metrics calculated before and after both harmonization approaches. Age correlated with unharmonized network metrics, but was more strongly correlated with network metrics based on both harmonization approaches. Parameter harmonization successfully controlled for scanner variability while preserving network topology and connectivity weights, indicating that harmonization of global network parameters based on unharmonized adjacency matrices may provide optimal results. The current work supports the use of ComBat for removing multiscanner effects on global network topology.

1. Introduction

Network neuroscience has become a popular approach to characterize brain structure in vivo in healthy and clinical populations (1–3). The structural connectome can be mapped using diffusion-weighted MRI (2, 4–6), a non-invasive technique that is sensitive to white matter microstructure (7).

Pediatric mild traumatic brain injury (mTBI) is a prevalent global public health concern (8–11) that is characterized by subtle and diffuse alterations in brain tissue [reviewed in (12, 13)]. The neurobiology of pediatric mTBI remains poorly understood [see (12)]. White matter microstructural alterations can occur after pediatric mTBI, and multiple studies have examined specific white matter tracts using diffusion tensor imaging [DTI; (13–15)]. Emerging evidence indicates that pediatric mTBI can alter global and regional brain networks (16–19). Thus, network neuroscience may be a potentially promising tool that could provide a robust characterization of network mechanisms involved in this important and highly prevalent neurological disorder.

Large, multisite neuroimaging studies of pediatric mTBI have become increasingly common to reduce false positive results from small samples, increase statistical power, and enhance reproducibility and generalizability of results (20, 21). For instance, the Advancing Concussion Assessment in Pediatrics (A-CAP) study (22) is the largest neuroimaging study of pediatric mTBI to date, with recruitment occurring at five children's hospitals across Canada including longitudinal MRI assessment using 6 different scanners. The A-CAP study has the potential to increase scientific and clinical knowledge about neurobiological outcomes in pediatric mTBI. However, using multiple MRI scanners introduces non-biological data variability due to different scanner systems, models, and sequence protocols, among other factors (23–26). Managing these non-biological sources of variability in multisite studies is increasingly important to generate accurate group-level inferences and enable detection of underlying biological phenomena (23).

ComBat is a widely used method for multisite (multiscanner) harmonization that originated from techniques used for genomics data (27). It is one of the most well-validated tools for multiscanner harmonization of structural neuroimaging data that makes no assumptions about the origin of scanner variation (26). ComBat implements a multivariate linear mixed effects regression with terms for biological variables and site to model the features of interest; the model parameters are estimated using an empirical Bayes approach. For diffusion tractography, ComBat has already demonstrated higher performance for multiscanner harmonization than other methods such as removal of artificial voxel effect by linear regression (RAVEL) and functional normalization of metrics (26). Unlike a general linear model approach that includes site or scanner as a fixed effect covariate, ComBat demonstrates better outlier robustness to account for small within-scanner sample sizes by borrowing information across features to shrink estimates toward a common mean (27, 28). The multiplicative scanner effects are also corrected by removing heteroscedasticity of model errors across scanners (29). Furthermore, ComBat preserves the variability contributed by true biological effects [e.g., sex and age; (26)]. However, no study has yet examined whether ComBat is suitable for graph theory metrics derived from the structural connectome based on DTI.

Unlike tractography, which yields a final value for each white matter tract, connectome analyses use weighted adjacency matrices to calculate network parameters. In tractography or region-of-interest analyses, multisite harmonization is performed on final metrics [e.g., average fractional anisotropy measures; (26, 30–32)]. However, graph theory analysis takes place after connections in a network have been mapped, mathematically represented as an adjacency matrix, and summarized by the computation of network parameters (33). Two distinct approaches to data harmonization are therefore possible in network analysis: (1) before the calculation of network parameters (i.e., matrix harmonization; harmonization at the level of connectivity weights), or (2) after calculation of network parameters (i.e., parameter harmonization). Identifying the optimal timing of data harmonization during data processing and analysis may influence the harmonization of multisite data, and hence has important implications for the accuracy of conclusions drawn from multisite connectivity studies.

To our knowledge, the performance of ComBat for multiscanner harmonization in studies of network topology and neurological disorders has not been evaluated. Therefore, the present study examined the application of two approaches to data harmonization across sites in a sample of DTI scans from children with mTBI or mild orthopedic injury (OI).

2. Methods

2.1. Study Design and Procedure

Data were drawn from the Advancing Concussion Assessment in Pediatrics (A-CAP) study (22), a multisite prospective, cohort study with longitudinal follow-up in children [Mean age (range) = 12.37 ± 2.34 years (8.00–16.99 years); 289 (59.7%) Male; see Table 1] with pediatric mTBI (n = 313) or mild orthopedic injury (OI; n = 171). Briefly, children were recruited within 48 h of injury from five children's hospitals across Canada (32), all of which were members of Pediatric Emergency Research Canada [PERC; (34)], and returned for three post-injury follow-up assessments: post-acute (targeted for 10 days post-injury), 3 and 6 months. Injuries and acute signs and symptoms were assessed during an initial emergency department visit that took place within 48 h post-injury.

TABLE 1

Table 1. Demographic information for the participants at each site/scanner.

The study was conducted with the approval of the research ethics board at each study site. All participants provided written informed assent and parents/guardians provided written informed consent (22). This study examined data from the MRI scans collected during the post-acute visit, as previously described (32, 35).

2.2. Diffusion MRI

Eligible participants completed a 3T MRI scan without sedation at the post-acute visit [see (35) for details]. In brief, thirty diffusion-weighted images with different diffusion gradient encoding directions were acquired at b = 900 s/mm², along with five images at b = 0 s/mm², with 2.2 mm isotropic resolution at all sites [General Electric: TR/TE = 6, 12 s/70, 90 ms; Siemens: 6.3, 7.8 s/55, 90 ms; (22)]. Data collected in Montreal was acquired using two different scanners, coded as Montreal 1 and Montreal 2, for a total of 6 scanners (“sites”): Calgary (General Electric), Edmonton (Siemens), Montreal 1 (General Electric), Montreal 2 (Siemens), Ottawa (Siemens), and Vancouver (General Electric).

2.2.1. Quality Assurance

Visual quality checks of all raw images were conducted to identify and exclude scans with structural abnormalities/incidental findings, scanner artifacts (e.g., warping), incomplete acquisition, or not collected using the standardized scan parameters (32). Data that passed the initial quality assessment were subsequently rated for motion by two trained analysts. Discrepancies were resolved through a third reviewer blind to initial ratings. Diffusion-weighted volumes with severe motion artifact were removed, and any scans with > 7 volumes with severe motion artifact were excluded from subsequent analysis (36).

2.2.2. Structural Connectome

Detailed image processing methodology has been previously described (18). Briefly, ExploreDTI (37) was used to preprocess diffusion images, calculate the diffusion tensor, conduct whole brain fiber tractography, and compute an adjacency matrix for each participant. Preprocessing included correction for signal drift (38), eddy currents, subject motion with rotation of the B-matrix (39), and susceptibility distortions (40). A deterministic streamline approach was used for whole brain fiber tractography (randomized seed points; seed and tractography FA threshold = 0.10; step size = 0.50 mm; angle threshold = 30°; step size = 0.5 mm; streamline length 50–500 mm). The resulting whole brain fiber tractography was extracted and used to compute an adjacency matrix for each participant.

The automated anatomical labeling [AAL-90, (41)] template was used to define 90 nodes in native (diffusion) space using functions from open-source software packages in MATLAB R2019a [see (18)]. Fully connected 90 × 90 adjacency matrices were constructed using the average FA of passing fibers among nodes in ExploreDTI for each participant and an absolute threshold of 0.10.

2.3. Global Network Metrics and Multisite Harmonization

The following global network parameters were calculated in MATLAB using the GRETNA software toolbox [http://www.nitrc.org/projects/gretna/; (42)]: global efficiency, global clustering coefficient, small worldness, modularity, and density. Network parameters were normalized against 1,000 randomly generated matrices.

Global network parameters were evaluated before harmonization and after two different harmonization approaches: matrix harmonization and parameter harmonization. The steps used for each approach are summarized in Figure 1. For both harmonization approaches, ComBat v1.0.5 (https://github.com/Jfortin1/ComBatHarmonization/tree/master/R) was conducted in R v3.6.3 [(43); https://www.R-project.org/] to harmonize the data for scanner differences. A covariate matrix with group (mTBI/OI), age at injury, and biological sex was included to preserve this variance:

$\begin{array}{l}mod<-model.matrix(~injury+age+sex)\end{array}$

FIGURE 1

Figure 1. Overall study procedure illustrating the data processing steps for the generation global network parameters (A) before harmonization, and the implementation of (B) matrix harmonization, and (C) parameter harmonization.

2.3.1. Matrix Harmonization

For matrix harmonization, weighted connectivity matrices were harmonized for multiple scanners and global network parameters were calculated using the harmonized connectivity matrix for each participant (see Figure 1). First, the lower diagonal values of each connectivity matrix were extracted to construct a dataframe of 4,005 columns corresponding to node connection pairs among the 90 defined brain regions (nodes), excluding self-connections [i.e., principal diagonal; (n(n-1))/2]. This was done because undirected adjacency matrices are diagonally symmetrical. ComBat was then used to harmonize those extracted values:

$\begin{array}{l}neuroCombat(dat=LowerDiagonal,batch=Site,mod=mod)\end{array}$

After the harmonization of the extracted connectivity weights, the adjusted square and symmetric weighted matrices were reconstructed for each participant and subsequently used to calculate global network parameters. During matrix harmonization, many of the connection weights that were 0 before harmonization (i.e., indicating that no connection existed between two nodes) were transformed to negative values. To correct for this transformation, an additional masking step was applied to reassign negative weights to 0 prior to graph analysis. Specifically, the masking step multiplied the binary connectivity matrix derived before harmonization with the harmonized weighted connectivity matrix for each participant (see Figure 1).

2.3.2. Parameter Harmonization

For parameter harmonization, the raw global network metrics (i.e., calculated before harmonization) were harmonized for multiple scanners using ComBat. Each parameter was harmonized in separate models because the distribution of each parameter is not necessarily related to the distribution of other parameters. The empirical bias estimation option was not applied (i.e., eb = FALSE) during parameter harmonization because each global network parameter was harmonized separately (i.e., number of features < n). For each global network metric, the following model was used to harmonize the data:

$\begin{array}{l}neuroCombat\\ (dat=Parameter,batch=Site,mod=mod,eb=FALSE)\end{array}$

2.4. Statistical Analysis

Statistical analyses were conducted using R v3.6.3. To evaluate the performance of each harmonization approach (i.e., adjacency matrix or network parameters), the effect of site was examined using separate one-way ANOVA models for each global network parameter. Non-significant scanner effects (p > 0.05) were interpreted as a successful removal of variability due to different scanners.

The proportion of significant (p < 0.05) post-hoc pairwise between-site comparisons was evaluated by calculating the number of significant uncorrected pairwise comparisons across scanners, divided by the total number of possible pairwise comparisons (i.e., n = 15). Correction for multiple comparisons was not applied for post-hoc t-tests followups, providing a more conservative evaluation of scanner effects.

The within-scanner consistency of the global network metrics before (unharmonized) and after each harmonization approach (matrix harmonization, parameter harmonization) was examined by calculating the intraclass correlation coefficient (ICC), with ICC < 0.50, 0.50 ≤ ICC < 0.75, 0.75 ≤ ICC < 0.90, and ICC ≥ 0.90 indicative of poor, moderate, good, and excellent consistently, respectively (44). Successful harmonization would reduce the effect of site while preserving the within-scanner variability for each parameter observed before harmonization.

To evaluate whether ComBat harmonization preserves biological variability, analysis of covariance (ANCOVA) was used to examine the effect of site, age at injury, sex, and group (mTBI, OI) on each network parameter. Significant effects involving age at injury were further examined using Pearson correlations, which were compared using a back-transformed average Fisher's Z procedure for dependent and overlapping correlations (45), as implemented using the cocor package (46). Overlapping correlations were used to conduct the following pairwise comparisons of age correlations: (1) matrix harmonization vs. unharmonized data, (2) parameter harmonization vs. unharmonized data, and (3) parameter harmonization vs. matrix harmonization.

Within-scanner age correlations on the unharmonized data were calculated to provide a reference value for the expected age correlation for each network parameter following harmonization. The reference value was calculated based on the means of within-site age correlations, weighted by the corresponding sample size of each scanner. Weighted means were calculated because sites with a greater number of participants may influence the correlation values to a greater extent than sites with smaller cohorts. Successful preservation of age-related biological variability across all scanners following harmonization would approximate the weighted mean of within-scanner age-correlations. Group differences between mTBI and OI were calculated using t-test.

3. Results

3.1. Presence of Site/Scanner Effects Before Harmonization

Before harmonization, all global network metrics differed by site (see Table 2 and Figures 2A, 3A). The largest site effect was observed for global efficiency [F₍₅₎ = 651.08, p < 0.001], with 14 of 15 (93%) significant between-site comparisons, followed by modularity [F₍₅₎ = 309.87, p < 0.001; 13 (86%) significant between-site comparisons], density [F₍₅₎ = 286.23, p < 0.001; 13 (86%) significant between-site comparisons], small worldness [F₍₅₎ = 182.93, p < 0.001; 12 (80%) significant between-site comparisons], and clustering coefficient [F₍₅₎ = 86.38, p < 0.001; 12 (80%) significant between-site comparisons].

TABLE 2

Table 2. Results summarizing the overall effect of site on global metrics before harmonization, after matrix harmonization and after parameter harmonization.

FIGURE 2

Figure 2. Violin plots illustrating the distribution of values across sites for global network parameters calculated (A) before harmonization, after (B) matrix, and (C) parameter harmonization.

FIGURE 3

Figure 3. Heatmaps illustrating pairwise between-site differences and t-values (lower diagonal) and within-site ICC values (principal diagonal) for the global network parameters calculated (A) before harmonization, after (B) matrix, and (C) parameter harmonization.

3.2. Matrix Harmonization

Main effects of site remained significant for global network metrics after matrix harmonization (Table 2 and Figures 2B, 3B). However, pairwise site differences were less pervasive after harmonization for global efficiency [F₍₅₎ = 3.88, p < 0.001; 4 (26%) significant between-site comparisons] and small worldness [F₍₅₎ = 170.69, p < 0.001; 11 (73%) significant between-site comparisons], but not for modularity [F₍₅₎ = 158.06, p < 0.001; 13 (86%) significant between-site comparisons], density [F₍₅₎ = 286.23, p < 0.001; 13 (86%) significant between-site comparisons], or clustering coefficient [F₍₅₎ = 295.44, p < 0.001; 15 (100%) significant between-site comparisons].

Within-site consistency of unharmonized and matrix harmonized metrics ranged from poor to excellent (Figure 3B). The highest consistency was observed for density [ICC = 1.00], which is the only parameter that measures the presence of connections but ignores their weights. Global efficiency and modularity also had excellent ICCs [Mean (range) = 0.96 (0.91, 0.99) for both], whereas consistency was good for small worldness [Mean (range) = 0.85 (0.57, 0.95)] and poor for clustering coefficient [Mean (range) = 0.21 (0, 0.62)].

3.3. Network Parameter Harmonization

Site was not significantly associated with global network metrics after parameter harmonization (Table 2 and Figures 2C, 3C). The within-site pre-post harmonization ICCs (Figure 3C) were consistently excellent for global efficiency [Mean (range) = 0.99 (0.99, 1)] and density [Mean (range) = 0.98 (0.97, 1)], and good to excellent for clustering coefficient [Mean (range) = 0.97 (0.86, 1)] and modularity [Mean (range) = 0.91 [0.85, 0.99)] and moderate to excellent for small worldness [Mean (range) = 0.86 (0.73, 0.99)].

3.4. Relationships Between Network Topology and Age Before and After Harmonization

Age significantly correlated with the following unharmonized network parameters (Figure 4A): global efficiency (r = 0.16, p < 0.001), clustering coefficient (r = 0.13, p < 0.003) and density (r = 0.12, p < 0.006). After matrix harmonization (Figure 4B), the correlation between age and global efficiency (r = 0.38, p < 0.001) was larger than before harmonization (z = 3.63, p < 0.001). After parameter harmonization, age-correlations increased (Figure 4C) for efficiency (r = 0.44, p < 0.001), clustering coefficient (r = 0.21, p < 0.001), and density (r = 0.27, p < 0.001), although only age correlations for efficiency (z = 4.71, p < 0.001) and density (z = 2.40, p < 0.016) were significantly larger compared to the unharmonized data.

FIGURE 4

Figure 4. Scatter plots illustrating the Pearson correlations between age at injury and each global network parameter calculated (A) before harmonization, after (B) matrix, and (C) parameter harmonization.

No significant correlations with age were observed for modularity and small worldness before harmonization, but modularity significantly correlated with age following matrix harmonization (r = −0.11, p =0.018), and both parameters showed significant age relationships following parameter harmonization (modularity: r = −0.15, p < 0.001; small worldness: r = −0.13, p < 0.003), although the coefficients were not significantly higher compared to the unharmonized data (p > 0.05). Clustering coefficient (z = 0.26, p < 0.008) and density (z = 2.40, p < 0.016) were significantly stronger following parameter as compared to matrix harmonization. Within-group age correlations are reported in Supplementary Table 1. No significant effects of sex or group were observed for any of the network parameters (see Supplementary Table 2).

4. Discussion

The popularity of large, representative datasets from collaborative, multisite research initiatives and structural connectomics has increased in recent years. Previous studies demonstrated that ComBat can control for site (scanner) differences while preserving biological variability (e.g., due to injury group, age, sex) for several MRI modalities (26, 29, 47, 48). This is the first study to validate the use of ComBat for the structural connectome. Here, ComBat was successfully used to harmonize structural connectivity data based on diffusion-weighted MRI across multiple scanners (“sites”). Parameter harmonization reduced the variability associated with different scanners to a greater extent than matrix harmonization, although both approaches reduced site differences in global network metrics. As expected, both harmonization approaches also preserved biological effects of age on network parameters. Moreover, expected age-related associations with global network parameters were stronger after applying parameter as opposed to matrix harmonization. Overall, the results extend the validity of using ComBat harmonization to network parameters derived using diffusion-weighted MRI.

Parameter harmonization showed superior performance for removing scanner effects compared to matrix harmonization. Furthermore, parameter harmonization is more computationally efficient. Matrix harmonization requires a series of steps that involve value translation. Specifically, connectivity weights were deconstructed from the matrices by extracting the lower (or upper) diagonal elements, organized in a high dimensional data frame for harmonization, and reconstructed back in square matrices following harmonization. In some instances, this approach transformed connection weights that were initially 0 (i.e., no connection exists between two nodes in unharmonized data) to negative values, requiring an additional step reassigning these values to zero before graph analysis. Thus, matrix harmonization preserved the location, but not the strength of connections among node pairs. In contrast, parameter harmonization requires only one step. This appears beneficial in preserving the true global properties of the network, as illustrated by the reduced variability of the within-site consistency between the parameter harmonized and unharmonized global network metrics (see Figure 3).

Before harmonization, global efficiency exhibited more robust site effects than other measures, such as clustering coefficient. Matrix harmonization reduced (e.g., global efficiency), introduced (e.g., clustering coefficient), or maintained (e.g., density) site effects compared to the unharmonized data. The variable performance of matrix harmonization across different metrics may indicate that properties of the network other than the pairwise connection strengths are affected by scanner. Except for density, global parameters included in the present analysis encode information about the topology (i.e., location) as well as the weights of connections among distinct brain regions. Density, which reflects the number of connections regardless of their strengths, did not demonstrate differences in the magnitude of site effects following matrix harmonization (see Figures 2, 3), indicating that site differences are present in topological properties of the network beyond the strengths of pairwise connections (e.g., the number or location of connections). In addition, clustering coefficient quantifies segregation across brain regions (i.e., nodes) by counting the occurrence of existing connections between groups of three nodes. Since matrix harmonization does not alter the location of connections, groups of connected nodes maintain their configuration before and after harmonization. Furthermore, the magnitude of scanner effects might differ slightly among connections, and matrix harmonization might differently impact the reciprocal connection strengths across groups of nodes (i.e., it targets pairs of nodes), potentially explaining the variable performance for removing site effects in the case of clustering coefficient (see Figures 3A,B). These topological properties may be better controlled by parameter harmonization, because global parameters already encode this information.

Correlations between age and network parameters were generally larger following both harmonization approaches, but were slightly more robust following parameter harmonization. One exception was the relationship between age and clustering coefficient, which weakened following matrix harmonization. This is in line with the other results, suggesting that matrix harmonization may be problematic for clustering coefficient. The detection of significant age effects following parameter harmonization, even in the absence of significant correlations for the unharmonized data, raises the question of whether additional variability was added to the data during harmonization that might have artificially boosted the relationship between age and network topology. Further analyses suggest this is not the case, because the age correlation following parameter harmonization were closer to the weighted means of within-scanner correlations before harmonization (see Supplementary Table 1). In addition, previous studies show relationships between age and global network topology in typical development (49–51) and children with TBI (19). This indicates that parameter harmonization may better preserve age-related biological variability compared to matrix harmonization and to the unharmonized data, although differences between the two harmonization approaches were small.

The children with mTBI and OI did not differ in any global network metrics before or after each harmonization approach. This was expected given that DTI and NODDI indices of white matter microstructure did not differ between groups previously in this sample (18, 32, 35), and other pediatric samples at similar time points (52). Another study compared a subset of this sample (children recruited at the Calgary site) to typically developing children and also did not find global or regional (nodal) network differences between mTBI and mild OI groups post-acutely, but did find an effect of injury more generally relative to typical development (18).

The current study did not address the effect of data harmonization applied prior to the generation of adjacency matrices, which is an additional possibility to account for the variability across different scanners [e.g., using methods described by (26, 53, 54)]. It has been suggested that connectome generation can be stable across scanners based on the derived network parameters (55). While future studies may consider this, data harmonization prior to connectome construction is increased in complexity, involving additional processing steps. These include warping the data into a common space and deconstructing brain images to build a voxel by participant data frame, which does not allow for the construction of adjacency matrices in native diffusion space. Following voxelwise harmonization, data would need to be reconstructed into subject-specific brain images (i.e., harmonized FA maps), which may impose substantial feasibility challenges due to the high computational complexity and number of additional transformations involved in this process.

There are some limitations to the current work. Weighted connectivity matrices were analyzed in this study; future multisite studies might examine whether differences in binary matrices relate differently to the effects of site. We did not assess the influence of different thresholds on harmonization. In addition, the current study used only one parcellation for the construction of adjacency matrices, and future studies might focus on whether other parcellations are similarly affected by site effects particularly when running matrix harmonization. While most methods use similar preprocessing steps, slight variations in these steps and how they are applied can impact calculated diffusion metrics, and thus may be important to explore in future studies. Data acquisition in the current study included single shell diffusion-weighted data. Multishell acquisition protocols may be differently affected by site effects, which might be addressed in future studies. Lastly, the current study used deterministic tractography, and future analyses might consider testing the effect of harmonization on networks derived using probabilistic tractography, as the two approaches have been shown to differ in terms of within- and between-scanner consistency (55).

5. Conclusions

The present paper validates the utility of ComBat harmonization in the context of graph theoretical analysis for structural connectivity derived from DTI. The harmonization of global parameters derived from unharmonized adjacency matrices provided superior performance as compared with the harmonization of connectivity weights for removing between- site differences, preserving the within-site variability and preserving age-related biological variability in the data.

Data Availability Statement

A dataset with deidentified participant data and a data dictionary will be made available upon reasonable request from any qualified investigator, subject to a signed data access agreement.

Ethics Statement

The current study was reviewed and approved by each site where data was collected. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

Author Contributions

AO, AW, AH, KY, and CL: conception and design of the study, acquisition and analysis of data, and drafted significant portion of the manuscript and figures. MB, CB, WC, QD, SF, and RZ: conception and design of the study. All authors contributed to the article and approved the submitted version.

Funding

Funding provided by the Canadian Institutes of Health Research (FDN143304), Ronald and Irene Ward Chair in Pediatric Brain Injury (KY); Canada Research Chair (CL, AH, and CB); Doctoral Fellowship from the IMT School for Advanced Studies Lucca (AO); Harley N. Hotchkiss-Samuel Weiss and Killam Postdoctoral Fellowship (AW); Alberta Children's Hospital Foundation Professorship in Child Health and Wellness (SF).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fneur.2022.850642/full#supplementary-material

References

1. Bassett DS, Sporns O. Network neuroscience. Nat Neurosci. (2017) 20:353–64. doi: 10.1038/nn.4502

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bullmore E, Sporns O. Complex brain networks: graph theoretical analysis of structural and functional systems. Nat Rev Neurosci. (2009) 10:186–98. doi: 10.1038/nrn2575

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Stam CJ. Modern network science of neurological disorders. Nat Rev Neurosci. (2014) 15:683–95. doi: 10.1038/nrn3801

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Gong G, He Y, Concha L, Lebel C, Gross DW, Evans AC, et al. Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cereb Cortex. (2009) 19:524–36. doi: 10.1093/cercor/bhn102

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Iturria-Medina Y, Canales-Rodríguez EJ, Melie-García L, Valdés-Hernández PA, Martínez-Montes E, Alemán-Gómez Y, et al. Characterizing brain anatomical connections using diffusion weighted MRI and graph theory. Neuroimage. (2007) 36:645–60. doi: 10.1016/j.neuroimage.2007.02.012

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Iturria-Medina Y, Sotero RC, Canales-Rodríguez EJ, Alemán-Gómez Y, Melie-García L. Studying the human brain anatomical network via diffusion-weighted MRI and Graph Theory. Neuroimage. (2008) 40:1064–76. doi: 10.1016/j.neuroimage.2007.10.060

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Lerch JP, van der Kouwe AJ, Raznahan A, Paus T, Johansen-Berg H, Miller KL, et al. Studying neuroanatomy using MRI. Nat Neurosci. (2017) 20:314–26. doi: 10.1038/nn.4501

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Cassidy JD, Carroll L, Peloso P, Borg J, Von Holst H, Holm L, et al. Incidence, risk factors and prevention of mild traumatic brain injury: results of the WHO Collaborating Centre Task Force on Mild Traumatic Brain Injury. J Rehabil Med. (2004) 36:28–60. doi: 10.1080/16501960410023732

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Thurman DJ. The epidemiology of traumatic brain injury in children and youths: a review of research since 1990. J Child Neurol. (2016) 31:20–7. doi: 10.1177/0883073814544363

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Centers for Disease Control and Prevention [CDC]. Report to Congress on Traumatic Brain Injury in the United States: Epidemiology and Rehabilitation. National Center for Injury Prevention and Control (2015). p. 1–72.

Google Scholar

11. Yeates KO. Traumatic brain injury. In: Yeates KO, Ris MD, Taylor HG, Pennington BF, editors. Pediatric Neuropsychology: Research, Theory, and Practice. New York, NY; London, UK: Guilford Press (2010). p. 112–46.

Google Scholar

12. Mayer AR, Kaushal M, Dodd AB, Hanlon FM, Shaff NA, Mannix R, et al. Advanced biomarkers of pediatric mild traumatic brain injury: progress and perils. Neurosci Biobehav Rev. (2018) 94:149–65. doi: 10.1016/j.neubiorev.2018.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Lindsey HM, Hodges CB, Greer KM, Wilde EA, Merkley TL. Diffusion-weighted imaging in mild traumatic brain injury: a systematic review of the literature. Neuropsychol Rev. (2021) 1–80. doi: 10.1007/s11065-021-09485-5

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Schmidt J, Hayward KS, Brown KE, Zwicker JG, Ponsford J, Van Donkelaar P, et al. Imaging in pediatric concussion: a systematic review. Pediatrics. (2018) 141:e20173406. doi: 10.1542/peds.2017-3406

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ware AL, Shukla A, Goodrich-Hunsaker NJ, Lebel C, Wilde EA, Abildskov TJ, et al. Post-acute white matter microstructure predicts post-acute and chronic post-concussive symptom severity following mild traumatic brain injury in children. NeuroImage. (2020) 25:102106. doi: 10.1016/j.nicl.2019.102106

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Imms P, Clemente A, Cook M, D'Souza W, Wilson PH, Jones DK, et al. The structural connectome in traumatic brain injury: a meta-analysis of graph metrics. Neurosci Biobehav Rev. (2019) 99:128–37. doi: 10.1016/j.neubiorev.2019.01.002

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Sharp DJ, Scott G, Leech R. Network dysfunction after traumatic brain injury. Nat Rev Neurol. (2014) 10:156–66. doi: 10.1038/nrneurol.2014.15

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Ware AL, Yeates KO, Geeraert B, Long X, Beauchamp MH, Craig W, et al. Structural connectome differences in pediatric mild traumatic brain and orthopedic injury. Hum Brain Mapp. (2021) 43:1032–46. doi: 10.1002/hbm.25705

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Watson CG, DeMaster D, Ewing-Cobbs L. Graph theory analysis of DTI tractography in children with traumatic injury. NeuroImage. (2019) 21:101673. doi: 10.1016/j.nicl.2019.101673

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Poldrack RA, Baker CI, Durnez J, Gorgolewski KJ, Matthews PM, Munafó MR, et al. Scanning the horizon: towards transparent and reproducible neuroimaging research. Nat Rev Neurosci. (2017) 18:115–26. doi: 10.1038/nrn.2016.167

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Maas AI, Menon DK, Adelson PD, Andelic N, Bell MJ, Belli A, et al. Traumatic brain injury: integrated approaches to improve prevention, clinical care, and research. Lancet Neurol. (2017) 16:987–1048. doi: 10.1016/S1474-4422(17)30371-X

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Yeates KO, Beauchamp M, Craig W, Doan Q, Zemek R, Bjornson BH, et al. Advancing concussion assessment in pediatrics (A-CAP): a prospective, concurrent cohort, longitudinal study of mild traumatic brain injury in children: study protocol. BMJ Open. (2017) 7:e017012. doi: 10.1136/bmjopen-2017-017012

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Pinto MS, Paolella R, Billiet T, Van Dyck P, Guns PJ, Jeurissen B, et al. Harmonization of brain diffusion MRI: concepts and methods. Front Neurosci. (2020) 14:396. doi: 10.3389/fnins.2020.00396

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhu T, Hu R, Qiu X, Taylor M, Tso Y, Yiannoutsos C, et al. Quantification of accuracy and precision of multi-center DTI measurements: a diffusion phantom and human brain study. Neuroimage. (2011) 56:1398–411. doi: 10.1016/j.neuroimage.2011.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhu AH, Moyer DC, Nir TM, Thompson PM, Jahanshad N. Challenges and opportunities in dMRI data harmonization. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. Shenzhen: Springer (2019). p. 157–72. doi: 10.1007/978-3-030-05831-9_13

CrossRef Full Text | Google Scholar

26. Fortin JP, Parker D, Tunç B, Watanabe T, Elliott MA, Ruparel K, et al. Harmonization of multi-site diffusion tensor imaging data. Neuroimage. (2017) 161:149–70. doi: 10.1016/j.neuroimage.2017.08.047

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. (2007) 8:118–27. doi: 10.1093/biostatistics/kxj037

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Beer JC, Tustison NJ, Cook PA, Davatzikos C, Sheline YI, Shinohara RT, et al. Longitudinal combat: a method for harmonizing longitudinal multi-scanner imaging data. Neuroimage. (2020) 220:117129. doi: 10.1016/j.neuroimage.2020.117129

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Fortin JP, Cullen N, Sheline YI, Taylor WD, Aselcioglu I, Cook PA, et al. Harmonization of cortical thickness measurements across scanners and sites. Neuroimage. (2018) 167:104–20. doi: 10.1016/j.neuroimage.2017.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Kurokawa R, Kamiya K, Koike S, Nakaya M, Uematsu A, Tanaka SC, et al. Cross-scanner reproducibility and harmonization of a diffusion MRI structural brain network: a traveling subject study of multi-b acquisition. NeuroImage. (2021) 2021:118675. doi: 10.1016/j.neuroimage.2021.118675

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zavaliangos-Petropulu A, Nir TM, Thomopoulos SI, Reid RI, Bernstein MA, Borowski B, et al. Diffusion MRI indices and their relation to cognitive impairment in brain aging: the updated multi-protocol approach in ADNI3. Front Neuroinformatics. (2019) 13:2. doi: 10.3389/fninf.2019.00002

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Ware AL, Yeates KO, Tang K, Shukla A, Onicas AI, Guo S, et al. Longitudinal white matter microstructural changes in pediatric mild traumatic brain injury: An A-CAP study. Hum Brain Mapp. (2022) 1–15. doi: 10.1002/hbm.25885. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Bullmore ET, Bassett DS. Brain graphs: graphical models of the human brain connectome. Annu Rev Clin Psychol. (2011) 7:113–40. doi: 10.1146/annurev-clinpsy-040510-143934

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Bialy L, Plint A, Zemek R, Johnson D, Klassen T, Osmond M, et al. Pediatric emergency research Canada: origins and evolution. Pediatr Emerg Care. (2018) 34:138–44. doi: 10.1097/PEC.0000000000001360

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Shukla A, Ware AL, Guo S, Goodyear B, Beauchamp MH, Zemek R, et al. Examining brain white matter after pediatric mild traumatic brain injury using neurite orientation dispersion and density imaging: an A-CAP study. NeuroImage. (2021) 2021:102887. doi: 10.1016/j.nicl.2021.102887

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ware AL, Shukla A, Guo S, Onicas A, Geeraert BL, Goodyear BG, et al. Participant factors that contribute to magnetic resonance imaging motion artifacts in children with mild traumatic brain injury or orthopedic injury. Brain Imaging Behav. (2021) 1–12. doi: 10.1007/s11682-021-00582-w

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Leemans A, Jeurissen B, Sijbers J, Jones D. ExploreDTI: a graphical toolbox for processing, analyzing, and visualizing diffusion MR data. In: Proceedings of the International Society for Magnetic Resonance in Medicine, Vol. 17. Hawaii (2009).

Google Scholar

38. Vos SB, Tax CM, Luijten PR, Ourselin S, Leemans A, Froeling M. The importance of correcting for signal drift in diffusion MRI. Magn Reson Med. (2017) 77:285–99. doi: 10.1002/mrm.26124

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Leemans A, Jones DK. The B-matrix must be rotated when correcting for subject motion in DTI data. Magn Reson Med. (2009) 61:1336–49. doi: 10.1002/mrm.21890

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Veraart J, Sijbers J, Sunaert S, Leemans A, Jeurissen B. Weighted linear least squares estimation of diffusion MRI parameters: strengths, limitations, and pitfalls. Neuroimage. (2013) 81:335–46. doi: 10.1016/j.neuroimage.2013.05.028

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Tzourio-Mazoyer N, Landeau B, Papathanassiou D, Crivello F, Etard O, Delcroix N, et al. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. Neuroimage. (2002) 15:273–89. doi: 10.1006/nimg.2001.0978

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Wang J, Wang X, Xia M, Liao X, Evans A, He Y. GRETNA: a graph theoretical network analysis toolbox for imaging connectomics. Front Hum Neurosci. (2015) 9:386. doi: 10.3389/fnhum.2015.00386

PubMed Abstract | CrossRef Full Text | Google Scholar

43. R Core Team (2020). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Available online at: https://www.R-project.org/

Google Scholar

44. Koo TK, Li MY. A guideline of selecting and reporting intraclass correlation coefficients for reliability research. J Chiropract Med. (2016) 15:155–63. doi: 10.1016/j.jcm.2016.02.012

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Hittner JB, May K, Silver NC. A Monte Carlo evaluation of tests for comparing dependent correlations. J Gen Psychol. (2003) 130:149–68. doi: 10.1080/00221300309601282

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Diedenhofen B, Musch J. COCOR: a comprehensive solution for the statistical comparison of correlations. PLoS ONE. (2015) 10:e0121945. doi: 10.1371/journal.pone.0121945

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Pomponio R, Erus G, Habes M, Doshi J, Srinivasan D, Mamourian E, et al. Harmonization of large MRI datasets for the analysis of brain imaging patterns throughout the lifespan. NeuroImage. (2020) 208:116450. doi: 10.1016/j.neuroimage.2019.116450

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Yu M, Linn KA, Cook PA, Phillips ML, McInnis M, Fava M, et al. Statistical harmonization corrects site effects in functional connectivity measurements from multi-site fMRI data. Hum Brain Mapp. (2018) 39:4213–27. doi: 10.1002/hbm.24241

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Chen Z, Liu M, Gross DW, Beaulieu C. Graph theoretical analysis of developmental patterns of the white matter network. Front Hum Neurosci. (2013) 7:716. doi: 10.3389/fnhum.2013.00716

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Koenis MM, Brouwer RM, van den Heuvel MP, Mandl RC, van Soelen IL, Kahn RS, et al. Development of the brain's structural network efficiency in early adolescence: a longitudinal DTI twin study. Hum Brain Mapp. (2015) 36:4938–53. doi: 10.1002/hbm.22988

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Lim S, Han CE, Uhlhaas PJ, Kaiser M. Preferential detachment during human brain development: age-and sex-specific structural connectivity in diffusion tensor imaging (DTI) data. Cereb Cortex. (2015) 25:1477–89. doi: 10.1093/cercor/bht333

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Wilde EA, Ware AL, Li X, Wu TC, McCauley SR, Barnes A, et al. Orthopedic injured versus uninjured comparison groups for neuroimaging research in mild traumatic brain injury. J Neurotrauma. (2019) 36:239–49. doi: 10.1089/neu.2017.5513

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Fortin JP, Sweeney EM, Muschelli J, Crainiceanu CM, Shinohara RT, Initiative ADN, et al. Removing inter-subject technical variability in magnetic resonance imaging studies. NeuroImage. (2016) 132:198–212. doi: 10.1016/j.neuroimage.2016.02.036

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The SVA package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. (2012) 28:882–3. doi: 10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Bonilha L, Gleichgerrcht E, Fridriksson J, Rorden C, Breedlove JL, Nesland T, et al. Reproducibility of the structural brain connectome derived from diffusion tensor imaging. PLoS ONE. (2015) 10:e0135247. doi: 10.1371/journal.pone.0135247

PubMed Abstract | CrossRef Full Text | Google Scholar