As part of the TOP-NT consortium project, there were four data acquisition sites: the University of California Los Angeles (UCLA), the University of Florida and Morehouse School of Medicine (UF/MSM), Georgetown University and Uniformed Services University (GU/USU), and John Hopkins University (JHU), and a data analytics site—University of California San Francisco (UCSF). All study protocols used were in compliance with the Public Health Service Policy on Humane Care and Use of Laboratory Animals and approved by each site-specific Animal Research Committee, UCLA Approval #2021–082; UF/MSM Approval # 202010145; GU/USU Approval APG-21-044; JHU Approval #RA22M370. A total of 186 adult male and female, Sprague Dawley rats were used for this study and were randomly assigned to each group (UCLA: 33 TBI, 17 sham, JHU: 35 TBI, 14 sham, GU/USU: 33 TBI, 16 sham, UF/MSM: 27 TBI, 11 sham). Diffusion imaging data were acquired at 3 and 30 days after either controlled cortical impact (CCI) injury or sham anesthetic controls in adult, male (Mean ± SD body weight 317 ± 43.7 g, 271.3 ± 20.3 g, 225 ± 14.8 g, 369 ± 16.1 g, at UCLA, GU/USU, JHU and UF/MSM, respectively) and female rats (body weight 254 ± 19.1 g, 239.3 ± 13.8 g, 183 ± 12.7 g, 256 ± 7.9 g, at UCLA, GU/USU, JHU and UF/MSM, respectively). A total of 343 scans were collected across the two post-injury times across all sites (25 datasets were missing). Prior to data acquisition, the injury and imaging protocols were harmonized across the sites using standard operating procedures (7, 19). Post-acquisition data analysis, including all pre- and post-processing data steps were also harmonized using a BASH script that wrapped multiple commands which was shared across the sites. The detailed methods below document the cross-site harmonized methods and any differences across sites are noted. All data were acquired with adherence to the ARRIVE guidelines (20) including group assignment randomization, blinding by use of unique animal IDs across site with semi-automated analysis routines, as well as conforming to the other eight guidelines on protocol, analysis and reporting standards.

Rats were housed in pairs before and after the surgery in standard cages with ad libitum access to water and food. Rats were maintained on a 12:12 h light/dark cycle in the vivarium. Estrus cycle was assessed for female rats by vaginal lavage/swab at the time of surgery to determine the stage of the estrus cycle of the rat on the injury day. Two different levels of CCI were induced at each site by varying the deformation depth below the dura of the impactor using similar methods across the sites (Table 1). Rats were anesthetized with 5% isoflurane vaporized in 1 L/min O2 (2% isoflurane was used for maintenance of the anesthesia with adjustments between 1.5 and 2%), and then were head-fixed in a stereotactic frame and their body temperature was maintained at 37°c using a temperature-controlled heating-pad. After prophylactic administration of analgesics according to the site-specific animal protocol, aseptic technique was used to make a midline, longitudinal skin incision, followed by a 6 mm craniectomy (5 mm UF/MSM) centered at +3.5 mm anteroposterior and 3 mm lateral to the midline over the left hemisphere using a dental drill cooled under intermittent saline. CCI injury was conducted using a 5 mm (4 mm UF/MSM) diameter, bevel-edge, metal impactor tip at an angle parallel to the dural surface and at impact depth below dura, speeds and dwell times ranging from: 1.3–2.3 mm, 3.5-6 m/s and 100-240 ms, using an electromagnetic impact system (Leica Biosystems, IL, United States). The craniectomy was closed with a non-toxic sealant (Kwik-Cast, WP instruments, United States) and the skin sutured.

Data were acquired on Bruker Biospec consoles (Bruker, Billerica, MA, United States) connected to magnets from a variety of vendors and with different field strengths (Tesla, T) across the sites (7 T -UCLA and GU/USU, 11.7 T – JHU and 11.1 T—UF/MSM) with gradients of a varying peak strength and rise time resulting in different slew rates. Radiofrequency coils were similar across three sites (birdcage volume transmit coil decoupled from a receive-only 4-channel phased-array surface coil) and a transceive quadrature coil at UF/MSM. Rats were imaged prone under isoflurane sedation (1–1.5%) vaporized in medical air (0.5 L/min) positioned in a cradle with head stabilization using a 3-point head fixation. The diffusion data were acquired after a number of other scans within the same imaging session as part of a larger protocol, taking around 60 min prior to diffusion imaging. The prior scans for functional and structural imaging were conducted under a mixture of continuous dexmedetomidine sedation (0.05 mg/kg/h, subcutaneously) and isoflurane (0.5%) prior to switching to isoflurane alone. Temperature was maintained at 37°C via homeothermic-controlled external air or water heating.

A 3-dimensional, single-shot, spin echo, echo planar imaging sequence was used to acquire diffusion-weighted images with directionally encoded, monopolar diffusion gradients applied along 88 different, non-colinear directions, with b values of 0, 1,000 and 3,000 s/mm² (n = 4, 44 and 44, respectively), and no averaging. Four left–right, phase-reversed, b = 0 images were also acquired. Repetition and echo times varied slightly according to the supported gradient hardware (Table 2). Similarly, the diffusion gradient amplitude (big delta) was varied across site to arrive at similar b values while maintaining similar diffusion times, the time between application of the diffusion-sensitizing gradients (Table 2). The acquisition data matrix was constrained to a 72x49x96 matrix size, in the 1st phase-encoding, 2nd phase-encoding and read-out directions, within a field of view (FOV) of 18×12.25x24mm, respectively in the dorsal-ventral, left–right and antero-posterior planes, resulting in an isotropic resolution of 250 μm. Outer-volume signal suppression using multiple, 5-10 mm saturation slices were applied to reduce signal from outside the brain.

Whole brain anatomical data were acquired using a 3D-multi-gradient echo sequence with a variable matrix size and field-of-view across site, resulting in the same isotropic resolution of 160 μm (Table 2). Eleven echoes were acquired from 2.8–51.8 ms with an echo spacing of 4.9 ms, a repetition time of 45-125 ms and a 13–20 flip angle (Table 2).

Bruker raw data were converted to NIFTI (21), brain extracted (22), denoised (23), and corrected for Gibbs ringing (24), after which phase distortions were unwrapped with TOPUP (25, 26) and Eddycorrect (27), all of which were implemented under MRTRIX (28). Data were then fit to the tensor to derive scalar images- fractional anisotropy (FA), axial, radial and mean diffusivity (AD, RD, MD, respectively). FA data were then used to derive an unbiased mean deformation template (MDT) at each site using successive registrations of rigid, affine and non-linear alignment under ANTS (29). The resulting affine transformations and warp fields were applied to all scalar data to align them to the local site MDT. Voxel-wise sham mean and standard deviation (SD) maps were computed for each tensor indices, which were then used together with each co-registered brain volume scalar to calculate regions of high and low indices relative to local sham data at a corrected z-score threshold based on a distribution-corrected z-score (30, 31) beginning with a desired z value of |z| > ±3.1 (p < 0.001). This DisCo-Z correction is required to adjust for group-wise bias originating for differences in sample size (30). The resulting corrected z values for thresholding the injury groups at each site at 03d/30d was: 3.85/3.75, 3.99/3.85, 4.19/4.19 and 5.08/5.08, respectively at UCLA, GU/USU, JHU, UF/MSM. For the sham group the adjusted z values at 03d/30d was: 2.67/2.72, 2.61/2.67, 2.52/2.52 and 2.38/2.38, respectively at UCLA, GU/USU, JHU, UF/MSM. The brain-wide volumes of tissue surviving these upper and lower z-score thresholds for each tensor scalar, herein designated Scalar_LOW and Scalar_HIGH volumes, were then quantified for each brain by counting the voxels that were greater or less than the z-score. For voxel-level harmonization, a multi-site MDT was calculated using the site-specific MDTs as input, and then each site-specific subject space data were transformed to the mean site MDT in a single step through application of the original warp and affine transformation together with those generated from site MDT to multi-site MDT space. Univariate harmonization statistics (e.g., corrected z-scores) were recalculated in the multi-site MDT space using the same method that had been applied at the site-specific level. Tensor-Based Deformation Analysis: Anatomical data were used to assess local tissue deformations to approximate tissue atrophy/compression and swelling/expansion as we have before (32, 33). Briefly, a mean deformation brain template was constructed from all data using ANTS (34) and the Jacobian Determinant was calculated from the resulting warp transformation fields. Sham Jacobian data were used to calculate voxel-based mean and standard deviation maps which were then used to calculate z statistic maps for all injured data. The volume of tissue surviving a threshold of z ± 3.01 (p < 0.01, uncorrected) was determined for each injured brain and used to assess the region of local tissue contraction and referred to herein as tissue atrophy.

The R package NeuroCombat_1.0.13 was used for harmonization.¹ The covariates included in the model were day post-injury, atrophy, group, and sex (Figure 1). Site was used as the batch parameter. Given prior successes in applying NeuroCombat to cross-sectional data (scans less than 180 days apart) (16, 35) and its ability to adjust for covariates (36), all scans across days and groups were combined. By using large amounts of data in the model (samples >20), we hoped to replicate prior successes correcting for site effects in site-to-site traveling subject data (36). Atrophy was used as a proxy for the level of injury rather than the surgeon-designated mild/moderate scale. This is because designated injury severity is still not well modelled. Using atrophy as a continuous variable also afforded greater sensitivity to differences compared to a binary choice of gross injury severity. Days post injury were defined as the number of days since anesthesia in shams, and this was treated as a continuous variable. Empirical Bayes was left enabled as per the default NeuroCombat settings. All statistical analyses were conducted in RStudio version 2023.6.0.421 (Posit Team 2023). Cohen’s d was calculated as a measure of effect size by dividing the difference between the mean injured and mean sham animal by the pooled standard deviation. The False Discovery Rate (FDR) correction was applied to correct for multiple comparisons of all statistics.

The number of brain voxels with FA greater or less than a statistical theory-derived distribution-corrected z-score (equivalent to p < 0.001 relative to site-specific sham data, Figure 2) are referred to as FA_High and FA_Low whole-brain volumes of pathology, respectively. We have operationally defined “whole-brain volume of pathology” based upon a voxel tensor scalar value being lower or greater than the z-score threshold, as before (31). It should be noted that variability among measurements in shams is not due to injury, but other likely sources of variance, such as between scanner effects and animal-to-animal differences. Relevant code to reproduce the analyses in this study can be found at the following link: https://github.com/ngharris/neurocombat_TBI.

Univariate Harmonization: The volume of brain regions significantly different from sham was derived from harmonization of univariate volume data derived from subject data co-registered to site-specific brain template space and is herein referred to as univariate-level harmonization (Figure 1). Ggplot2 was used for plotting violin plots for univariate data (37). Four-site data was visualized with residual plots of values predicted by a regression model generated from the harmonized data across each site to demonstrate the presence of unequal variance among sites post-harmonization. As described above, NeuroCombat harmonization used atrophy, day post-injury, group and sex as covariates. Site was used as the batch parameter.

Since the data obtained from brain volumes used for univariate harmonization are count data derived from the number of voxels above or below a z-threshold, we also implemented Combat-seq. Combat-seq was originally designed for RNA sequencing count data, but can be natively applied to neuroimaging count data (18). We used sex, atrophy, and day post injury as covariates and treatment as the group variable. Site was used as the batch parameter. Bayes shrinkage remained disabled. Due to non-normality in the data, univariate metrics were compared using Kruskal Wallis tests to test for effect of group and site, and were FDR corrected for multiple comparisons. The estimation stats package was used to visualize effect size (38). p-values less than 0.05 were considered statistically significant for univariate statistics. Group differences were also maintained after harmonizing AD, MD, and RD data.

Previous work (43) has shown that if the proportion of outliers is higher than 5% of the data, or if extreme outliers exist, this can severely distort the ability of NeuroCombat to remove site effects. In the current data, we found that the variance for both high and low FA values was significantly increased within Site 4 sham data as compared to all other sites for univariate data before harmonization (p < 0.001). This increased variance within site 4 data was maintained even after harmonization when compared to the other three sites (p < 0.001). Consistent with this data, the residuals of FA values for sham rats predicted by a regression model generated with harmonized sham data showed unequal distributions with the largest spread in site 4 data. One potential cause of this heteroscedasticity could be attributed to the presence of outliers, which could lead to ineffective harmonization if all sites were included (Figure 3). However, log transformation of these data revealed similar results with a greater number of univariate outliers in site 4 compared to all and significantly increased variance. Removal of site 4 data resulted in an absence of site effects in the harmonized data (p > 0.05), further substantiating these findings. As a result, all further univariate analyses were conducted with data from only three sites in order to more accurately determine the utility of using NeuroCombat to harmonize preclinical imaging data from multiple sites. Additional analyses using data from all four sites can be found in Supplementary Figures S16–S18.

An important consideration when investigating differences in multi-site data is whether to pool shams from all sites, or whether to treat site-specific sham data independently for the calculation of univariate volume data. Pooled shams were used to generate the z-scored FA_Low and FA_High values for each rodent prior to harmonization. To gauge how effectively differences between sham and TBI rats were detected when considering pooled sham data versus sham data specific to each site, we examined FA_Low and FA_High distributions across the sites. Compared to unharmonized data, harmonization with pooled shams led to a 32% increase in the number of injured rats across all sites, where FA_Low volumes of pathology were more than two standard deviations larger than the mean sham animal across all sites (n = 21 vs. 30 animals at day 3 and n = 19 vs. 23 at day 30, respectively). The number of injured rats with FA_High volumes above the mean after harmonization were reduced or remained similar to unharmonized data (n = 9 vs. 8 at day 3, respectively and no change at day 30).

Harmonization with site-specific sham data afforded a lower improvement compared to the prior reported pooled sham harmonization; the number of injured animals with FA_Low volumes above sham mean data after harmonization increased by only 10% (n = 33 vs. 32 at day 3 and n = 52 vs. 45 at day 30, respectively). This aligns with previous suggestions that “mega-analyses,” in which pooling is done prior to statistical analyses, results in an improved ability to detect injury compared to meta-analyses that do not pool data (12). However, the number of injured rats with significantly different FA_High volumes after harmonization increased by 17.6% (n = 11 vs. 9 at day 3 and n = 9 vs. 8 at day 30, respectively). While univariate-level harmonization led to mixed results overall, the greatest improvement occurred in injured rat FA_Low volumes using pooled sham data. As a result, since NeuroCombat is intended for use on data pooled from multiple sites, pooled shams were used for all further univariate analyses. Pooled shams were also used for voxel-level analyses for consistency.

Previous reports have demonstrated that removal of site-specific effects after harmonization would indicate that any remaining group differences could accurately be attributed to a true difference in pathology, thus demonstrating successful harmonization (12). Within the original, unharmonized, univariate data, there were significant sham vs. TBI group effects for FA_Low volume day 3 data (Figure 4A, p < 0.001, FDR correction, Kruskal-Wallis test) and significant site effects for FA_High volume at both time points (p < 0.001, Figures 4B,D). After univariate-level harmonization with Combat-seq, the group difference in FA_Low volume originally found at day 3 was retained, and a group difference at day 30 was introduced (p < 0.001, Figures 4A,C). Finally, univariate-level harmonization failed to correct for the site effects observed in FA_High volumes at day 3 and 30 post-injury before harmonization (Figures 4B,D) where sham data from Site 1 had substantially higher volumes of pathology and were more variable than the other two sites (p < 0.001, Kruskal-Wallis test). This result was similarly found in the other measures, where group differences were maintained post-harmonization, but a significant effect of site was no longer found (Supplementary Figures S1, S6, S11).

Despite these unexpected effects, we calculated the effect size for each site in order to quantify any improvement in the ability to discriminate between experimental groups due to harmonization. Using the difference in univariate FA volume for injured versus pooled sham data at each timepoint, effect size changed in different directions for day 3 FA_Low data after harmonization at sites 1, 2 and 3 by −4.3, 1.35 and 2.1-fold, respectively. This was significantly different from zero (p < 0.05) for Sites 2 and 3 compared to pre-harmonization, when it was not previously different (Figure 5A). Similar but smaller effects were found for day 30 data at all sites (Figure 5B). Despite the introduction of group differences in the FA_High data at day 3 through harmonization, the magnitude of the effect size decreased for site 3 (p < 0.05 pre-harmonization, p > 0.05 after harmonization, Figure 5C) and it remained similar for other sites (p > 0.05).

Unexpectedly, neither Combat-seq or NeuroCombat were able to completely resolve the site effects in univariate FA_Low data on day 3, however, Combat-seq was able to uncover strong group effects at day 30 (Figures 4A,C; data not shown for NeuroCombat). NeuroCombat was also less sensitive for detecting FA_Low injury compared to Combat-seq, but similar for FA_HIGH (data not shown for NeuroCombat). This occurred at both the day 3 and day 30 timepoints for FA_Low, where Combat-Seq showed a more significant effect of group compared to NeuroCombat (p < <0.001 versus p < 0.05, respectively, Kruskal Wallis test). Group differences after harmonizing AD, MD, and RD data using Combat-seq were maintained with both NeuroCombat and Combat-seq.

By visual inspection of the group-level incidence maps of volumes of pathology described by FA_Low and FA_High, we found there were marked differences at site 4 compared to Sites 1–3 at both 3 and 30 days after injury (Figure 2), in agreement with the prior outlier volume analysis above. We confirmed this quantitatively for FA values; of the four sites analyzed, 10.5% of the FA values from Site 4 were found to lie more than 1.5 interquartile ranges outside the interquartile range of FA values when compared to the other three sites at each voxel location. The outliers among Sites 1, 2 and 3 were smaller, containing only 4.6, 1.7 and 0.9% outliers, respectively, which was more than 2-fold less than Site 4. Additionally, Site 4 FA population data were found to have an excess kurtosis of 0.497, compared to −0.774, 0.159 and 0.274 in Sites 1, 2, and 3, respectively indicating a quite different distribution of FA volumes. These discrepancies were similar to the findings at the univariate volume level described in the prior sections. Harmonization of whole brain data resulted in large decreases in the standard deviation of Site 4 FA data and corresponding increases in variance in Sites 2 and 3, as well as to a lesser extent in Site 1 (Figure 6A). Since the presence of an outlier site severely biased the performance of NeuroCombat, further analysis of whole brain data was done with only Sites 1, 2 and 3.

Bland–Altman plots were used to quantify the effect of harmonization across the whole brain between sites at the voxel level across all subjects (Figures 6B,C). These graphs plot the mean value of the same voxel location from each pair of sites against their difference, and they were generated for each pair of sites before and after harmonization. In the original, unharmonized data, large differences were observed for all combinations of sites indicating site-related differences. After NeuroCombat harmonization, the difference between voxel measures between pairs of sites approached zero, indicating that FA values were more similar at the same brain location across sites after harmonization. This was supported by a Kolmogorov–Smirnov analysis, which showed significantly different site effects before harmonization (p < 0.001) and a resolution of these site effects after harmonization on average between each combination of sites (p > 0.05, Figures 6B,C).

Crucially, all group level effects shown prior to harmonization through univariate volume analysis were maintained using voxel-level harmonization: group difference in FA_Low at 3- and 30-days post injury were maintained (p < 0.01), as well as site-specific differences in FA_High at 3- and 30-days post injury (p < 0.01) (Figure 4).

Since the application of NeuroCombat at the voxel level resulted in decreased cross-site differences, supporting prior work at the same site across different field strengths (17), we next investigated the statistical effect of harmonization within regions of primary injury as well as in less injured regions across the brain. Correct interpretation of statistical results requires consideration not only of statistical significance but effect size and statistical power. Improvements in power, the probability of detecting a result when it is present, can help to detect significant differences in datasets with small effect sizes. Given the simple statistical design of the current study, increases in sample size from combining data in this study should generally be expected to increase statistical power at a geometric rate through reduction in standard error. However, power is also affected by effect size and unexplained variances, among other factors. Technical/hardware variation and site related effects when measuring the injury present themselves as ideal sources of variance for NeuroCombat to capture. We aimed to understand how cross-site harmonization was able to resolve this variance, as well as how harmonization affected effect size and power across the brain. We hypothesized that harmonization would increase the probability of correctly detecting group-level differences, particularly within specific brain regions that may exhibit high variability at the single site level through increases in power and effect size.

To accomplish this, whole brain images of change in power and effect size due to harmonization were generated by voxel-based subtraction between the original, non-harmonized and the harmonized FA data for each site (Figure 7). The average effect from harmonization calculated over all three sites was an increase in power around the primary injury area of the ipsilateral corpus callosum but decreases in most other regions (Figure 7A). At the site level at day 3, power increased within the ipsilateral white matter and gray matter surrounding the primary injury area in all three site data, with the greatest improvement measured in Site 1 (Figure 7B). There were, however, decreases in power as a result of harmonization, globally within the grey matter across all sites and at both times post-injury. At day 30, the greatest increase in power was observed in data from Site 1, similar to day 3, while there were power decreases in the striatum and ventral areas of the brain in data from Sites 2 and 3. These same remote changes also generalized to the AD, MD and RD (Supplementary Figures S4, S9, S14).

Calculation of the corresponding voxel-wise change in effect size revealed an ipsilateral increase in the magnitude in the site-averaged data (Figure 8A) consistent with the prior power changes (Figure 7A). However, while power changes were decreased by harmonization throughout much of the brain, effect size increased over the majority of the brain remote from the primary injury site, with the exception of decreases within multiple clustered voxels that occurred at both post-injury times. Across site-averaged data, these trends persisted in AD, MD and RD (Supplementary Figures S4, S9, S14), but the drivers varied. For FA, the underlying data driving these site-wide effects were a decrease in effect size at site 1, but an increase at Sites 2 and 3 (Figure 8B). The opposite trend was seen in MD and RD, where Site 1 showed increases in effect size in the injury area, while Sites 2 and 3 exhibited large decreases around the injury area (Supplementary Figures S9, S14). AD showed a different trend, with Site 2 driving much of the decreases and Sites 1 and 3 showing increased effect size at the injury area (Supplementary Figure S4). The heterogeneity in both effect size and power across the sites outside of the primary injury area may be explained by the effect of NeuroCombat as it attempts to adjust for the variability introduced by site effects, ultimately converging at a common spatial pattern.

We sought to determine whether the improved power and effect size was underpinned by the detection of pathologically low or high FA when compared to the sham group. FA voxel overlap maps showing the proportion of rats with FA values significantly different from sham after harmonization versus before, revealed an increased number, and overlap of rats after harmonization when compared to the original data (Figure 9). These increases were especially prevalent within, and adjacent to, the primary, ipsilateral injury area. All differences detected were related to decreases in FA; no data was underpinned by increased areas of FA_High when compared to shams. These results indicate that cross-site data harmonization can improve the reliability of group differences in the areas surrounding the primary injury area at the cost of lower power in regions more remote from the primary injury site. For other scalar measures, results were markedly different. Generally, the proportion of rats significantly different from mean sham was lower. In the AD dataset, most areas of difference were found in the ipsilateral cortex, frontal pole and cerebellum. The only location to produce an increase in the number of rats was a small portion of the corpus callosum. In other areas, harmonization decreased the number of voxels with significant differences, suggesting an initially weak effect at both day 3 and 30 (Supplementary Figure S5). Similarly, the MD data at day 3 showed large portions of the cortex with a small number of rats being significantly different from sham, but this area decreased after harmonization. At day 30, the number of voxels increased around pre-existing areas of difference in the cerebellum and parts of the ipsilateral cortex (Supplementary Figure S10). In the RD dataset, areas of pre-existing cortical differences decreased after harmonization at day 3 and day 30, except for minor increases in coverage in the contralateral cerebellum (Supplementary Figure S15).

A significant challenge to NeuroCombat harmonization is the presence of outliers. These findings are consistent with prior reports (43) that the presence of outliers can severely bias the ability of NeuroCombat to remove site effects. That prior work identified outliers as measurements outside 3 interquartile ranges (IQRs) of Q1 or Q3, and determined that if more than 5% of values at a particular site were outliers, then the ability for NeuroCombat to harmonize data was hindered, and this was associated with reduced variability within the outlier site data (43). Site 4 in the current study exhibited these characteristics with an outlier proportion of >10% when outliers were measured as 1.5 IQRs outside of Q1 and Q3. Difference in scanner field strength is unlikely to be the source of variation since another site that operated at the same field strength (Site 3) showed no such variation. One major difference in hardware between Site 4 and all other sites was the use of different radiofrequency coil hardware that resulted in a reduced brain coverage that may have accounted for the increased variation (Figure 2). A second major difference was the lower number of rats included from Site 4 in both experimental groups, making the statistical correction to the z-scored volume data larger and thus less sensitive to the detection of abnormal diffusion scalar-related values. This suggests that a more conservative approach during pre-acquisition harmonization is necessary when selecting if site data can be grouped together for harmonization. This may relate to both MRI hardware and injury severity level.

Use of pooled sham data from across all three sites, rather than a site-specific approach, produced a greater percentage of TBI rats where the volume of pathology indicated by FA_Low or FA_High was greater than in sham rats. This was counter to expectations, where our initial hypothesis was that site-specific sham data would encode site-specific technical, scanner-related effects that would be negated through within-site identification of TBI volumes of pathology, leading to more uniform cross-site data. This premise was based on the current lack of available preclinical phantoms that could be used for correcting for site-specific effects, leaving us to make use of local sham data to fill this gap. However, despite the greater success of using pooled shams, it was only the use of site-specific shams that retained group differences in univariate-level harmonized FA_Low volumes. Although Combat-seq was unable to completely resolve site effects in some univariate data, it did maintain or enhance the strong group differences. Because a similar group effect was found through voxel harmonization, it may indicate that Combat-seq is able to improve injury detection within univariate data. This may suggest that it is more suitable compared to NeuroCombat for harmonizing univariate, count-distributed neuroimaging data. This is the first study that the authors are aware of that has applied Combat-seq to univariate statistics derived from neuroimaging data. Combat-seq was previously only applied to RNA sequencing data. Because the voxel-level harmonization data did possess site effects, harmonization was able to preserve true biological variation, which prevented the loss of signal. Despite this finding, voxel-level harmonization maintained all group differences present within the original data. It is possible that the lower sample size at the single-site level could lead to effects that may either overestimate or underestimate true effects, where pooling data can increase sensitivity to true effects (44). It could also be that the whole-brain data used to generate voxel-level harmonization metrics was continuous, which enabled more successful harmonization with NeuroCombat. It therefore remains crucial to study the effect of how sham data are used when determining the effects of cross-site harmonization.

Consistent with prior literature (12), voxel-level harmonization effectively removed site effects in whole-brain data. As shown in the Bland–Altman plots (Figures 6B,C), across all combinations of sites and mean FA value of voxel pairs, the difference between each pair of voxels was reduced after harmonization. This reduction in variability across sites after harmonization corresponds with results seen in both the preclinical and clinical spaces (12, 17, 45). These changes suggest that harmonization improves the ability to detect group differences in areas of pre-existing biological variation. In areas where no significant differences between sham and injury were initially found, power decreased. Univariate measures of pathology showed similar results, with mild effect size increases between injured and sham groups after harmonization. This could explain some of the detrimental effects of harmonization, which has been shown to obscure biological variation in data without site effects, such as in regions unaffected by injury or site (16). The current finding that the proportion of rats which were significantly different from shams did not increase substantially outside of brain regions where group differences were detected before harmonization suggests that NeuroCombat does not introduce erroneous biological variation (14, 15). Voxel-level harmonization appears to only be able to improve detection of group differences across data where group differences already exist, rather than revealing differences obscured by noise or systematic error.

In areas of the brain where group differences were already present, NeuroCombat successfully improved delineation of injury among FA datasets, but not MD, RD or AD. Improvements in power and effect size were seen in areas of the ipsilateral white matter and cortex, which generally corresponds to areas of existing group differences in the FA data. This was further seen by improvements in the proportion of rats that exhibited group differences at these voxels at day 3 and day 30. The use of NeuroCombat on data with site effects may help to improve reproducibility in this respect, as improved power, effect size and more frequent detection of injury can help to reduce erroneous conclusions. An increase in power indicates a benefit for pooling data across sites, allowing for smaller samples from each of multiple sites to support observed effects. The decrease in power observed in non-injured areas could have been related to a lack of biological effect in those areas prior to harmonization, where harmonization increased power in areas of existing difference at the expense of whole-brain differences. An alternative reason could be related to a decrease in inaccurate effect sizes within non-injured regions, which would cause power to decrease in those areas and reflect areas of true biological difference. Similar findings have been reported in clinical FA data: harmonization with NeuroCombat was able to increase the number of voxels associated with age in a clinical multi-site and multi-scanner study of DTI measures’ relationship with age (12), as well as improve power in the detection of differences among individuals with schizophrenia in T1 derived metrics of cortical thickness, surface area and subcortical volumes (5). The present study builds upon these initial findings by validating NeuroCombat’s efficacy in preclinical data, as well as by demonstrating that these improvements in the detection of group differences are confined within areas of existing biological variation.

There was a relative lack of success by harmonization to improve injury detection among AD, MD and RD datasets compared to FA, despite varying degrees of power and effect size improvements. Reasons for this are not clear, but it may relate to site differences in measurement of absolute diffusivities compared to a more equitable comparison of FA which is a normalized ratio of these diffusivities, and thus less likely to differ in absolute values across site. Further examination into whether normalization of unbounded metrics helps to improve the efficacy of NeuroCombat may be warranted to better understand this phenomenon. Other methods of ComBat that have been used in imaging data include generalized additive models (gamCombat) which accounts for non-linear covariates, longCombat which can be used to detect changes over time, and gaussian mixture model-combat that removes variation due to unidentified covariates (16, 46). NeuroCombat has been shown to be preferable when conducting cross-sectional data harmonization compared to longCombat (16). While in the present study we were focused on group-level comparisons at each timepoint, the use of longCombat may well be useful to capture differences in injury trajectories across site in future work.

Inherent in the preprocessing of image data are the potential errors associated with the requirement to conduct spatial transformations to ensure correct registration to a common 3D space. Multiple transformations could have impacted the quality and precision of the resulting measures of pathology. The heterogeneous nature of injury may also play a role in the decrease in power in non-significantly different brain areas. Due to potential minor differences in delivery of trauma, brain areas of some rats may have been more injured compared to similar areas in other rats. This may lead to less consistent measurements and interfere with NeuroCombat’s abilities to resolve site effects.

Univariate level volume analysis was conducted after identifying regions of Scalar_Low and Scalar_High through a statistical-based z-corrected threshold (30, 31). However, there is no agreed upon statistical correction required for the analysis of voxel-based data. It is possible that the strict threshold used here (p < 0.001) reduced the ability to accurately judge pathologic tissue within injured versus sham rats, as indicated by the low overlap between rats in site 4 (Figure 2). More work is required to determine the impact of statistical correction threshold on determining group differences after TBI.

The biological basis for FA_HIGH regions after TBI is uncertain, although it has been found before in the same spatially discrete patterns from FA_LOW regions (3). Noteworthy here though is that the benefit provided by multi-site harmonization was afforded by FA_LOW, and from which the biological underpinnings is far more certain.