Fusion in diffusion MRI for improved fibre orientation estimation: An application to the 3T and 7T data of the Human Connectome Project

Determining the acquisition parameters in diffusion magnetic resonance imaging (dMRI) is governed by a series of trade-offs. Images of lower resolution have less spatial specificity but higher signal to noise ratio (SNR). At the same time higher angular contrast, important for resolving complex fibre patterns, also yields lower SNR. Considering these trade-offs, the Human Connectome Project (HCP) acquires high quality dMRI data for the same subjects at different field strengths (3T and 7T), which are publically released. Due to differences in the signal behavior and in the underlying scanner hardware, the HCP 3T and 7T data have complementary features in k- and q-space. The 3T dMRI has higher angular contrast and resolution, while the 7T dMRI has higher spatial resolution. Given the availability of these datasets, we explore the idea of fusing them together with the aim of combining their benefits. We extend a previously proposed data-fusion framework and apply it to integrate both datasets from the same subject into a single joint analysis. We use a generative model for performing parametric spherical deconvolution and estimate fibre orientations by simultaneously using data acquired under different protocols. We illustrate unique features from each dataset and how they are retained after fusion. We further show that this allows us to complement benefits and improve brain connectivity analysis compared to analyzing each of the datasets individually.

. Phantom used for numerical simulations. Two crossing bundles were used to simulate white matter (ground-truth orientations are shown color-coded, red for left/right, green for up/down). These bordered with areas that simulated cerebrospinal fluid (CSF) and areas that simulated grey matter (GM).
First, the relaxation-free, noise-free signal at each WM voxel was simulated using a mixture of signals obtained from diffusion tensors. Specifically, in non-crossing regions an isotropic tensor and an anisotropic axially symmetric tensor were used (with mean diffusivity d WM for both and eigenvalues (λ 1 , λ 2 ) for the latter as shown in Table 1). The volume fraction for the anisotropic tensor was f sum =0.6. In crossing regions, two axially symmetric tensors were used with volume fractions 0.3 each, oriented along the crossing orientations. For the GM and CSF regions isotropic tensors were used (with mean diffusivities d GM and d CSF respectively). To mimic small differences in the diffusivities of each tissue type caused by the different acquisition protocols (differences in pulse duration δ, diffusion time Δ, maximum gradient strength, echo time, for instance (Qin et al., 2009)), we computed the DTI mean diffusivity of each tissue type, using the b=1000 s/mm 2 shell from the 3T and 7T in-vivo data and used these values in the simulations (see Table 1).
Relaxation effects were then added and the final signal S was: with S nr the simulated signal without any relaxation effects. TE and TR, as well as the b values and gradient directions used for the simulations were the ones of the actual 7T and 3T acquisitions (see Data Acquisition section). The relaxation times at different field strengths were obtained from the literature (Table 1). For the 7T, the T 2 for WM and GM were obtained from (Yacoub et al., 2003) and for CSF from (Wyss et al., 2013), while T 1 values were found in (Rooney et al., 2007). For the 3T, the T 2 for WM and GM were obtained from (Stanisz et al., 2005) and for CSF from (Neelavalli et al., 2009), while T 1 values were found in (Wansapura et al., 1999) for GM and WM and in (Lin et al., 2001) for CSF. Finally, to reflect the different spatial resolutions at each field strength (7T are high resolution (HR), 3T are low resolution (LR)), the version of the phantom corresponding to the 3T parameters was downsampled by a factor of two in each dimension. This represents a scenario of having a ratio of 2 between the voxel size at HR and LR, as was done with the real data. To simulate Rician noise, zero-mean Gaussian signal was added in quadrature. For a given signal to noise ratio (SNR) at HR, a scaled SNR level was used for the LR data. The scaling was found using the following formula for a 2D EPI sequence (Haacke et al., 1999): with FS being the field strength, V the voxel volume, N PE the number of phase encoding steps, BW Pixel the bandwidth per pixel in the readout direction and R the in-plane parallel imaging acceleration factor. Using the parameter values given in Table 1 and V LR =8V HR , we obtained a scaling of SNR LR =3.7585 SNR HR .

Effect of ignoring relaxation in the data fusion model
The RubiX model ignores tissue relaxation effects and their difference at different field strengths. Using the simulator we explored the influence of this simplification on the estimates. We compared the accuracy of RubiX estimates in two scenarios, one where relaxation is not considered for signal generation and another where relaxation effects are considered as described above (dark blue vs orange bar plots in Figure S2. A case without relaxation but higher SNR is also shown as a baseline (light blue bar plots). Ignore for now the green bar plots, as they are explained in the next section). As expected, the performance of the model is best (smallest errors) for the data that match the model assumptions best, i.e. the data simulated without any relaxation effects. However, in the more realistic scenario of considering relaxation, RubiX also performs very well. The orientation estimation works almost identically between the two cases (with the exception of a small error increase of ~7% for orientations at WM/CSF borders, Fig. S2A). The main performance penalties are observed in other parameters (the diffusivity, the anisotropic volume fractions and the baseline signal intensity, Fig. S2, panels C-E) and at tissue interfaces (mostly at WM/CSF borders). The relative errors increase by 12%, 15%, 78% at WM/GM interfaces and 28%, 25%, 113% at WM/CSF interfaces for the diffusivity, anisotropic volume fraction and S 0 intensity respectively. However, in all cases, the absolute errors remain within reasonable limits (smaller than 5-10%). Therefore, the simplification made in the RubiX model regarding relaxation effects seems to be a reasonable one, particularly if the estimation of orientations is of interest.

Which spatial model?
The spatial model used in RubiX is given by Eq. (6). This models the signal attenuation at the lower resolution as the attenuation of the sum of the signals at the high resolution voxels that intersect a low resolution voxel. Alternatively, another model uses a weighted sum of HR attenuations to represent the attenuation at the lower resolution . We evaluated the robustness of these two alternatives, when relaxation effects (and differences due to field strengths) are taken into account. Figure  S2 orange and green bar plots represent errors using the two models. The performance of the two models is very similar within WM, but their difference is prominent at the WM/CSF interfaces, where the latter model is much more sensitive and prone to error (as before the orientation estimates seem to be less sensitive, with mean errors less than 10 o even at the interfaces, but large errors can occur in the other parameters). We therefore used the spatial model of Eq. (6) in our data analysis. Figure S2. Errors in the estimates from the RubiX model on simulated data under different scenarios: when relaxation effects are ignored (dark blue), when relaxation effects are considered (orange), when relaxation effects are considered and an alternative spatial model is used (green). The data were simulated using SNR HR =10, except a baseline case where relaxation effects are ignored and SNR HR =50 (light blue). Mean and standard deviation of errors across indicated regions are shown in each case. A) Orientation error (in case of crossings, the mean error is shown). B) 95% cones of uncertainty for the estimated orientations (in case of crossings the mean is shown). C) Error in total anisotropic volume fraction, D) Error in Mean diffusivity. E) Errors in estimated baseline signal intensity S 0 .

Effect of mis-registrations
The success of the data fusion relies on an accurate alignment of the datasets to be fused. Deviations from a good alignment can introduce errors to the results, due to mixing of signals that correspond to different tissue locations; as they would do for traditional voxel-wise analysis, where diffusion-weighted volumes need to be aligned in the presence of distortions and head motion. We assessed the sensitivity of our framework to misalignments of the LR and HR datasets. With our current tools used for distortion correction and alignment (Glasser et al., 2013, Andersson andSotiropoulos, 2016), misalignments in the order of 0.25 voxels are expected as a worst-case scenario (see Fig. 4 in (Graham et al., 2016)). We simulated two scenarios with misalignments along different directions ( Figure S3). The misalignments were in the order of 0.25 LR voxels (or equivalently 0.5 HR voxels). Compared to a perfect alignment, the biggest absolute errors occur for the volume fractions and the diffusivities at the tissue interfaces. Orientation errors also increase, but even for the worst misalignment (and at SNR=10), the mean error remains less than 7 o . Figure S3. Errors in the estimates from the RubiX model on simulated data under different scenarios: when the LR and HR data are perfectly registered (orange), when the HR data are displaced half a voxel size along direction Y (grey) and when the HR data are displaced half a voxel size diagonally (along directions X & Y) (dark grey). Relaxation effects were considered. The data were simulated using SNR HR =10. Mean and standard deviation of errors across indicated regions are shown in each case. A) Orientation error (in case of crossings, the mean error is shown). B) 95% cones of uncertainty for the estimated orientations (in case of crossings the mean is shown). C) Error in total anisotropic volume fraction, D) Error in Mean diffusivity. E) Errors in estimated baseline signal intensity S 0 .

Simulations Summary
The main findings from our simulations are: a) Ignoring relaxation effects and acquisition timings in the generative model is a simplification. However, this simplification influences more profoundly the estimates of diffusivity and volume fractions at the WM-CSF interface ( Figure S2). Fibre orientations estimates are considerably less affected.
b) The spatial model proposed here (Eq. 6) behaves well at tissue interfaces and across different acquisition protocols ( Figure S2), where integration of signals with potentially large differences may occur.
c) The proposed framework relies on a very good alignment of the fused datasets. But even when this is not the case, the model still performs reasonably well ( Figure S3).

In-vivo Results
Similar to Figure 4 in the main paper, we present data and model predictions for inferior axial locations (where artefacts and inhomogeneities are worse, see frontal regions in Figure S4) and for a different subject with a larger brain size (which again inherently increases inhomogeneities particularly at high field, see Figure S5 and Figure  11). Even if the difference between model predictions and data is slightly higher in regions with more inhomogeneities, the trend is always the same with the model providing good predictions of both HR and LR datasets at all b values.  Top row shows the actual measurements (the 3T data have been downsampled, as used in the estimation). Middle row shows the predictions. To obtain model predictions, the mode of the estimated posterior distribution was used for each model parameter. Bottom row shows the difference of the prediction minus the measurement expressed as a percentage of the measurement. The pial surface is shown as a gray line. (HCP Subject ID: 109123, larger brain size than 102311 and therefore inherently larger inhomogeneities at high field). Supplementing Figure 7 in the main paper, we provide Figure S6 to show that certain exquisite features in the high-resolution data cannot be recovered by simple interpolation. We interpolated the 3T preprocessed data to the 7T spatial resolution (i.e. 1.05mm isotropic) using spline interpolation. We then estimated fibre orientations on the interpolated data using the voxelwise model and assessed how perpendicular or parallel to the white matter/gray matter boundary surface these orientations were (3T_spline). For reference we include the maps for all the other cases shown in Figure 7. As shown, even if the spline interpolation changes slightly the amount of orientations being perpendicular to the surface (e.g. area highlighted by the white arrow), it does not recover the features seen in the 7T data (noticeably more orientations perpendicular to the surface), which are preserved with RubiX data fusion. Figure S6. Angle between surface tangent and fibre orientations at the inflated WM/GM boundary surface. Orientations have been estimated as before using voxel-wise deconvolution on the 3T alone (estimates at 1.25mm isotropic resolution), voxel-wise deconvolution on a splineinterpolated version of the 3T data to the 7T resolution (estimates at 1.05mm isotropic resolution), voxel-wise deconvolution on the 7T alone (estimates at 1.05mm isotropic resolution), RubiX deconvolution on both 3T and 7T (estimates at 1.05mm isotropic resolution). For every surface vertex, the maximum dot product between fibre orientations (with volume fraction f>5%) at this location and the surface normal is computed. This is then converted to the colour-code angle shown on the inflated surface. To aid visualisation of these qualitative illustrations, Gaussian smoothing on the surface with 1mm FWHM was performed for all cases. (HCP Subject ID: 158035)