In this work, we propose a novel method to harmonize dMRI data acquired from multiple sites and scanners, for a joint analysis of the data for increased statistical power.

In our proposed pipeline, the inputs are dMRI images acquired at multiple sites. We set one of these sites as our reference site and the rest as our target sites. As output, we provide harmonized dMRI signals of the target sites, which are mapped to the reference site by removing inter-site variabilities but preserving intra-site variability. To account for regional-variations in the data, we use Freesurfer parcellated label maps, which also provide region-wise correspondences. The current work assumes that the acquired data has the same b-value and spatial resolution, even though the number of gradient directions could be slightly different. In our pipeline, we start by computing spherical harmonic (SH) coefficients of the dMRI signal S=[s₁...s_G]^T acquired along G unique gradient directions. In the SH basis, the signal S can be written as S=Σ_i,jY_ijC_ij, where Y_ij is a SH basis function of order i and phase j and C_ij are the corresponding SH coefficients. It is well-known that the l₂ norm of the SH coefficients at each order forms a set of rotation invariant (RISH) features⁵: |C_i|²=Σ_jC_ij².

Given the Freesurfer label maps for each subject, at each site, we compute the expected value of the region-wise RISH features as the average of the voxel-wise RISH features. As shown in Figure 1, these features vary significantly across sites as well as for different regions. Considering this fact, in the next step of our pipeline, we introduce a region based mapping Π such that average RISH features computed over all subjects in the target site match those in the reference site, i.e. E_t(Π(|C_i|²))=E_r(|C_i|²) . A simple Π, which satisfies this equality is Π(|C_i|²)=|C_i|²+E_r-E_t. Given this region-wise map of the RISH features, we need to find a voxel-wise mapping π to estimate harmonized SH coefficients:

Π(|C_i|²)=Σ_j [π(C_ij)]² =Σ_j|C_ij|² +E_r-E_t

We set π as a scaling function:

π(C_ij)=[(Π(|C_i|²))/(|C_i|²)] C_ij

The computed voxel-wise scales for a subject is shown in Figure 2. Applying the above mapping over all the regions of the brain, we can re-compute the signal at each voxel as ζ=Σ_iΣ_jπ(C_ij)Y_ij and project it back to the reference site.

We validated our method on diffusion data acquired from seven different sites (including two GE, three Philips, and two Siemens scanners) on a group of age-matched healthy subjects. Scanner details and subject numbers for each site are reported in Table 1. Since all the data was from matched healthy subjects, we do not expect to see biological differences between the group of subjects from each site. Therefore, we hypothesize that the RISH or standard diffusion measures at the reference and the target sites are statistically not different. To validate our hypothesis, we used t-test to compute p-values for RISH features and standard single tensor based diffusion measures between the reference and the target sites; these p-values were computed before and after applying our harmonization method. Computed p-values for generalized FA (GFA) are shown in Table 2. It can be seen that all p-values-after harmonization are larger than 0.05, which indicates that we could not reject the null hypothesis (similar results are observed for all the RISH features and FA). We also computed the whole brain within site coefficient of variation (CV) in GFA, which did not change much after our harmonization. For independent validation, we also show results using tract-based spatial statistics before and after harmonization. In Figure 3, it can be seen that after data harmonization, all scanner specific group differences in the white matter were removed. Our experimental results demonstrate that, for nearly identical acquisition protocol across sites, scanner-specific differences can be accurately removed using the proposed method.

To the best of our knowledge, this is a first method that has addressed the issue of harmonizing dMRI data acquired from multiple sites⁶. Our proposed method is independent of compartmental modeling (e.g., tensor, intra/extra cellular compartments, etc.), subject-dependent, and region specific. It maintains inter-subject variability, removes scanner specific differences, and makes it feasible to do joint analysis to significantly increase the sample size and statistical power of neuroimaging studies.