Reducing Susceptibility Distortion Related Image Blurring in Diffusion MRI EPI Data

Diffusion magnetic resonance imaging (MRI) is an increasingly popular technique in basic and clinical neuroscience. One promising application is to combine diffusion MRI with myelin maps from complementary MRI techniques such as multi-parameter mapping (MPM) to produce g-ratio maps that represent the relative myelination of axons and predict their conduction velocity. Statistical Parametric Mapping (SPM) can process both diffusion data and MPMs, making SPM the only widely accessible software that contains all the processing steps required to perform group analyses of g-ratio data in a common space. However, limitations have been identified in its method for reducing susceptibility-related distortion in diffusion data. More generally, susceptibility-related image distortion is often corrected by combining reverse phase-encoded images (blip-up and blip-down) using the arithmetic mean (AM), however, this can lead to blurred images. In this study we sought to (1) improve the susceptibility-related distortion correction for diffusion MRI data in SPM; (2) deploy an alternative approach to the AM to reduce image blurring in diffusion MRI data when combining blip-up and blip-down EPI data after susceptibility-related distortion correction; and (3) assess the benefits of these changes for g-ratio mapping. We found that the new processing pipeline, called consecutive Hyperelastic Susceptibility Artefact Correction (HySCO) improved distortion correction when compared to the standard approach in the ACID toolbox for SPM. Moreover, using a weighted average (WA) method to combine the distortion corrected data from each phase-encoding polarity achieved greater overlap of diffusion and more anatomically faithful structural white matter probability maps derived from minimally distorted multi-parameter maps as compared to the AM. Third, we showed that the consecutive HySCO WA performed better than the AM method when combined with multi-parameter maps to perform g-ratio mapping. These improvements mean that researchers can conveniently access a wide range of diffusion-related analysis methods within one framework because they are now available within the open-source ACID toolbox as part of SPM, which can be easily combined with other SPM toolboxes, such as the hMRI toolbox, to facilitate computation of myelin biomarkers that are necessary for g-ratio mapping.


INTRODUCTION
Diffusion magnetic resonance imaging (MRI) is an increasingly popular method in neuroscience and clinical research. Diffusion MRI permits visualisation of the mobility of water molecules within brain tissue, and provides measures of microstructural integrity and anatomical connectivity (Conturo et al., 1999;Shimony et al., 1999;Le Bihan et al., 2001;Kaden et al., 2007). Knowledge of how different brain regions communicate is essential for understanding brain anatomy and its relationship with cognitive processes (Jbabdi et al., 2015;Hodgetts et al., 2017;Assaf et al., 2019;Movahedian Attar et al., 2020), and diffusion imaging in the clinic can be used for various purposes including the characterisation of strokes and tumours (Neumann-Haefelin et al., 1999;Bergui et al., 2001).
However, diffusion MRI data are generally acquired with echo planar imaging (EPI; Turner and Le Bihan, 1990) -a fast MRI acquisition technique that is prone to spatial distortions related to off-resonance effects (Jezzard and Balaban, 1995). It is, therefore, important to correct spatial distortions in diffusion data to provide as accurate a representation as possible of the true anatomy. Fortunately, various methods to address spatial distortions in diffusion MRI are available. Field inhomogeneities due to magnetic susceptibility variations can be measured with additional MRI sequences for field mapping (Jezzard and Balaban, 1995;Reber et al., 1998;Hutton et al., 2002) or estimated from the distorted images themselves using physical models (Chang and Fitzpatrick, 1992). A popular and efficient approach uses two EPI images with reversed phase-encoding directionsalso known as "blip-up" and "blip-down" images ( Figure 1A)to estimate the field map and correct the susceptibility-related distortions (Andersson et al., 2003;Weiskopf et al., 2005;Holland et al., 2010;Ruthotto et al., 2012;Breman et al., 2020; other susceptibly distortion correction methods are also available, e.g., Irfanoglu et al., 2015;Schilling et al., 2020). The issue remains, however, that even after distortion correction is applied, the corrected diffusion MRI images can suffer from additional image blurring artefacts because the compression of voxels leads to a loss of information. As illustrated in Figure 1, blip-up and blipdown images are often combined by calculating the arithmetic mean (AM) of the two images to provide a single diffusion image ( Figure 1B, top).
Avoiding residual spatial distortions and image blurring is particularly important when combining diffusion MRI with other quantitative MRI contrasts that can reveal additional information that would otherwise not be available (Mohammadi and Callaghan, 2018). For example, multi-parameter mapping (MPM; Weiskopf et al., 2013) is an emerging neuroimaging technique that provides a comprehensive approach to acquiring multiple quantitative biomarkers, namely, the longitudinal relaxation rate R 1 , proton density PD, effective transverse relaxation rate R * 2 , and magnetisation transfer saturation MT sat (Tabelow et al., 2019). These quantitative parameters contain specific and complementary information regarding the microstructural properties of brain tissue, such as the extent of myelination, iron and water concentration (e.g., Weiskopf et al., 2015;Kirilina et al., 2020). Combining diffusion MRI with MPM thus allows for the examination of microstructural information related to the neuronal conduction velocity in the form of the axonal g-ratio (Mohammadi et al., 2015a;Mohammadi and Callaghan, 2020). Furthermore, by combining diffusion MRI and MPM structural maps, it becomes possible to perform multivariate analyses with biomarkers sensitive to complementary tissue features (Draganski et al., 2011;Chowdhury et al., 2013). While various methods using reverse phase-encoding are available to efficiently reduce susceptibility-related distortions (for example, see Gu and Eklund, 2019), these implementations often focus only on diffusion MRI data. As such, to combine diffusion MRI with other quantitative MRI data, such as MPM for g-ratio mapping, the use of multiple toolboxes implemented in different software packages, is currently required.
FSL, together with TBSS for group analyses, is probably the most widely used package for diffusion MRI (Smith et al., 2004), with TOPUP for distortion correction (Andersson et al., 2003) and eddy for motion and eddy current correction . To examine the g-ratio, however, the diffusion data have to be combined with a quantitative MRI biomarker for the myelin compartment, requiring laboratory-specific in-house software that is not widely available (Stikov et al., 2015;Cercignani et al., 2017;Mancini et al., 2018). In another analysis package, SPM, 1 the ACID toolbox 2 can be used for distortion correction using Hyperelastic Susceptibility Artefact Correction (HySCO2; Ruthotto et al., 2012Ruthotto et al., , 2013Macdonald and Ruthotto, 2018) and motion and eddy current correction using ECMOCO (Mohammadi et al., 2010(Mohammadi et al., , 2015b. Moreover, the hMRI toolbox, also in SPM 3 (Tabelow et al., 2019), provides an easy-touse framework to generate quantitative structural maps from MPM data. SPM, therefore, has the advantage of being able to process both diffusion and MPM data, making SPM the only widely accessible software that contains all the processing steps required to perform group analyses of g-ratio data in a common space. However, limitations have been identified in the ability of the HySCO2 methodology to reduce susceptibility-related distortion (Gu and Eklund, 2019), highlighting a need for improvement. Moreover, the HySCO2 method uses the AM to combine the oppositely phaseencoded data after correction, meaning that blurring artefacts are typically present.
Here, we sought to improve the pipeline for the HySCO2based distortion correction, an approach we call consecutive HySCO, and to employ a different method to combine the corrected blip-up and blip-down images to reduce image blurring by using a weighted average (WA). We achieved this by downweighting regions that were squeezed (e.g., in the blip-up images) while up-weighting the same (stretched) regions in the opposite (blip-down) image. To determine regions that were stretched and squeezed, we used the Jacobian of the HySCO2 estimated field map. We compared the performance of the proposed WA approach to the standard AM combination used by HySCO2 in 1 www.fil.ion.ucl.ac.uk/spm/ 2 www.diffusiontools.com 3 www.hmri.info FIGURE 1 | Illustration (modified from Jones and Mahadevan, 2013) of the origin of blurring artefacts after combining distortion corrected diffusion data, and how blurring can be avoided. (A) Susceptibility-related distortions in two EPI images acquired with reverse phase-encoding directions (i.e., "blip-up" and blip-down" images) can lead to localised squeezing and stretching. Importantly, where one image is squeezed (red box, top) the reverse phase-encoded image (green box, bottom) is stretched. (B) After distortion correction and resampling, the two images are typically combined into an arithmetic mean (AM, red-green-dashed box, top). Here, we propose an alternative combination, the weighted average (WA, green box, bottom), where only the stretched boxes are combined, which retains the spatial information. Therefore, spatial information is lost in the AM combination leading to blurring artefacts whereas it is retained in the WA combination.
SPM. Where appropriate, we also benchmarked these changes against the widely available default pipeline in FSL (as used in the human connectome project; Glasser et al., 2013). Finally, given that the overall goal was to have all the processing steps necessary to perform g-ratio analysis within one (SPM) package, we also assessed whether there was any benefit of using the WA corrected diffusion data over the AM corrected data for g-ratio mapping.

Participants
Eighteen participants took part in the study. They were aged between 20 and 40 years old, reported no psychological, psychiatric, neurological or behavioural health conditions. The mean age of the sample was 30.0 years (SD = 6.07) and included 9 females and 9 males. Participants were reimbursed £10 per hour for taking part which was paid at study completion. All participants gave written informed consent and the study was approved by the University College London Research Ethics Committee.

Diffusion MRI Data Acquisition
Diffusion-weighted images were collected on Siemens Magnetom 3T TIM Trio systems using 32 channel head coils. The protocol used was the Multi-Band Accelerated EPI Pulse Sequence developed by the Centre for Magnetic Resonance Research at the University of Minnesota (R012ac, R013a on VB17 4 ; Feinberg et al., 2010;Xu et al., 2013). Acquisition parameters were: resolution = 1.7 mm isotropic; FOV = 220 mm × 220 mm × 138 mm; 60 directions with 6 interleaved b0 images for a total of 12 b0 images in each of the blip-up and blip-down phase encoding directions, echo time (TE) = 112 ms, repetition time (TR) = 4.84 s, with a multiband acceleration factor of 3. The sequence was performed 4 timestwice with b-values of 1,000 s/mm 2 and twice with b-values of 2,500 s/mm 2 . The first acquisition of each set of b-values was performed with phase-encoding in the anterior to posterior direction (blip-up), the second in the posterior to anterior direction (blip-down). The total acquisition time was 22 minutes.

Structural MRI Data Acquisition
Whole brain structural maps of magnetisation transfer (MT) saturation, at an isotropic resolution of 800 µm, were derived from an MPM quantitative imaging protocol Callaghan et al., 2015Callaghan et al., , 2019. This protocol consisted of the acquisition of three multi-echo gradient-echo acquisitions with either PD, T1 or MT weighting. Each acquisition had a TR of 25 ms. PD weighting was achieved with an excitation flip angle of 6 0 , which was increased to 21 0 to achieve T1 weighting. MT weighting was achieved through the application of a Gaussian RF pulse 2 kHz off resonance with 4 ms duration and a nominal flip angle of 220 0 . This acquisition had an excitation flip angle of 6 0 . The field of view was 256 mm head-foot, 224 mm AP, and 179 mm right-left. The multiple gradient echoes per contrast were acquired with alternating readout gradient polarity at eight equidistant echo times ranging from 2.34 to 18.44 ms in steps of 2.30 ms using a readout bandwidth of 488 Hz/pixel. Only six echoes were acquired for the MT weighted volume to facilitate the off-resonance pre-saturation pulse within the TR. To accelerate the data acquisition, partially parallel imaging using the GRAPPA algorithm was employed in each phase-encoded direction (anterior-posterior and right-left) with forty integrated reference lines and a speed up factor of two. Calibration data were also acquired at the outset of each session to correct for inhomogeneities in the RF transmit field (Lutti et al., 2010, Lutti et al., 2012. The total acquisition time was 27 min.

Diffusion MRI Distortion Correction Pipelines
Standard HySCO pipeline. The diffusion MRI data were first processed using the standard pre-processing pipeline available in the ACID toolbox (see text footnote 2) within SPM12 (see text footnote 1). Each diffusion image (regardless of b-value or blip-up / blip-down acquisition) was separately corrected first for motion and eddy current artefacts (Mohammadi et al., 2010), and then for susceptibility-related distortion artefacts using the HySCO2 module (Ruthotto et al., 2012(Ruthotto et al., , 2013Macdonald and Ruthotto, 2018). HySCO is a tool for the hyperelastic susceptibility correction of diffusion data that corrects geometric deformations and intensity modulations in EPI images caused by susceptibility-related field-inhomogeneities. To this end, HySCO takes a tailored variational image registration problem and incorporates a physical distortion model that is aimed at minimising the distance of two oppositely distorted images subject to invertibility constraints.
The new consecutive HySCO pipeline. As with the standard HySCO pipeline, each diffusion image (regardless of b-value or blip-up / blip-down acquisition) was first separately corrected for motion and eddy current artefacts and then for susceptibilityrelated distortion artefacts using the HySCO2 module. Next, the "mean_b0" of the blip-up data was co-registered to the "mean_b0" of the blip-down dataset using a rigid-body transformation (spm_coreg) and the estimated transformation was applied to the full blip-down dataset. Tensor fitting using all the b-values  was then performed separately on each of the distortion corrected blip-up and blipdown datasets to estimate Fractional Anisotropy (FA) maps. These FA images were then brain-masked using ACID. Finally, HySCO2 was repeated but now using the distortion corrected and brain-masked FA maps as input instead of b0 images; the second HySCO2 field map being consecutively applied to the "pre-corrected" diffusion MRI data (Figure 2). FSL pipeline. While our main interest was in comparing the two HySCO approaches, the diffusion MRI data were also processed based on the widely-used human connectome project pipeline (Glasser et al., 2013) in FSL 6.0.1 (Smith et al., 2004) to provide external benchmarking. The blip-up and blip-down data were separately corrected for distortions using TOPUP (Andersson et al., 2003), and this was then used as the input into eddy, with the distortion, motion and eddy-current correction performed together  with the replace outliers option enabled .

Diffusion MRI Phase-Encoding Combination Methods
For both the consecutive HySCO and FSL methods, the distortion-corrected data sets with opposite phase-encoding polarity were combined using the ACID toolbox within SPM12.
First, the consecutive HySCO blip-up and blip-down distortion corrected data were combined by calculating the AM of each blip-up and blip-down image pair.
Second, the consecutive HySCO blip-up and blip-down distortion corrected data were combined using a WA. This aimed to minimise information loss due to susceptibility distortion blurring induced by local spatial compression. To maximise the effective spatial resolution when combining the two images acquired with opposite phase-encoding directions, a region that was stretched in one of the images, increasing effective resolution, was up-weighted whereas the opposite image was down-weighted using a WA: being a Fermi function with two parameters (x 0 , k 0 ), tuning the sensitivity of the weighting of the Jacobian in a given voxel. Both parameters were heuristically optimised. Note that the Jacobian of the field map estimated in the initial HySCO step (Step 1, Figure 2) was used for weighting. The Jacobians were defined as follows: with b H being the HySCO field in units of displacement in the phase-encoding direction in mm and p being the phase-encoding direction -here the anterior-posterior direction. Finally, to provide a benchmark, the FSL blip-up and blipdown distortion corrected data were combined by calculating the AM of each blip-up and blip-down image pair.

Structural MRI Data Pre-processing
The structural MRI data were processed for each participant using the hMRI toolbox (Tabelow et al., 2019) within SPM12. The default toolbox configuration settings were used, with the exception that correction for imperfect spoiling was additionally enabled (Corbin and Callaghan, 2021). The output MT saturation map, used as the anatomically faithful structural image, quantified the degree of saturation of the steady state signal induced by the application of the off-resonance prepulse, having accounted for spatially varying T 1 times and RF field inhomogeneity (Helms et al., 2008;Weiskopf et al., 2013).
Each participant's MT saturation map was segmented into white matter probability maps using the unified segmentation approach (Ashburner and Friston, 2005), but with no bias field FIGURE 2 | The new consecutive HySCO pipeline. In step 1, as per the standard HySCO pipeline, the HySCO module is applied to the original b0 images, estimating the HySCO field (HySCO b0 ). In step 2, tensor fitting is performed on the distortion corrected blip-up and blip-down data (b0 b0 ) to estimate Fractional Anisotropy maps (FA b0FA ), which are then used in a second round of HySCO for estimating the residual susceptibility related distortions HySCO FA (b0FA). The HySCO field map is the displacement in the phase-encoding direction in mm.
correction (since the MT saturation map does not suffer from any bias field modulation) and using the tissue probability maps developed by Lorio et al. (2016).

Comparison of the Diffusion MRI Distortion Correction Methods
Analyses to investigate the differences between the diffusion MRI distortion correction methods were performed in native space. For each participant, two FA maps were estimated for each distortion correction method (standard HySCO, consecutive HySCO, FSL); one FA map from the distortion corrected blipup data and the other from the distortion corrected blip-down data. FA maps were estimated using the Diffusion Kurtosis Fit (DKI), utilising all the b-values, in the ACID toolbox. Each FA map was then co-registered to an individual's MT saturation white matter probability map using a rigid-body transformation. The difference between the individual blip-up and blip-down FA maps was calculated, restricted to a white-matter mask. The white-matter mask was generated by thresholding the MT saturation white matter probability maps above 0.9. From the resulting difference image, the root mean square over all voxels in the white matter mask was calculated for each individual for each distortion correction method.
Two-tailed paired t-tests, calculated in SPSS v25 and thresholded at p < 0.05, were then performed to test for significant differences between the root mean square for each of the distortion correction methods across participants. Effect sizes are reported as Hedge's g av , which is Cohen's d calculated for repeated measures and corrected for the positive bias caused by using sample estimates, using the resources provided by Lakens (2013). Paired t-tests were used instead of a repeated measures ANOVA to treat each distortion correction method independently.

Comparison of the Diffusion MRI Phase-Encoding Combination Methods
Analyses to compare the effects of the phase-encoding combination methods were also performed in native space. For each participant, three pairs of FA and b0 maps were estimated, one pair for each phase-encoding combination method (consecutive HySCO AM, consecutive HySCO WA, FSL). The FA and b0 maps were estimated using DKI fit in the ACID toolbox.
Multichannel segmentation, using the unified segmentation approach (Ashburner and Friston, 2005), but with the tissue probability maps developed by Lorio et al. (2016), was performed on each FA and b0 map pair to define the white matter tissue segments. These diffusion-derived white matter probability maps were then co-registered to each individual's MT saturation white matter probability map using a rigid-body transformation. This allowed for comparison between the white matter probability map defined from the distortion corrected diffusion images, and the white matter probability map from the MT saturation structural map. The percentage overlap between the diffusion white matter probability map and the structural white matter probability map was then calculated for each phase-encoding combination method.
As well as investigating the overlap of the diffusion and structural white matter probability maps as a whole, we also focused our analyses on regions that suffered from greater distortion. To identify these areas, for each participant, a map of the Jacobians was calculated from the first HySCO iteration. The co-registration transformations determined above were applied to each participant's Jacobian map to co-register the Jacobian and MT saturation structural images. The percentage overlap of the diffusion and structural white matter probability maps was calculated in regions of higher distortions. The masks for the aforementioned regions were determined by thresholding the Jacobian map at values higher than 10, 20, and 30% deviation from unity.
Two-tailed paired t-tests, calculated in SPSS v25 and thresholded at p < 0.05, were performed to test for significant differences between the percentage overlaps for each of the phase-encoding methods (consecutive HySCO AM, consecutive HySCO WA, FSL), at each level of distortion (whole white matter map, 10%, distortion, 20% distortion, 30% distortion). As with the distortion correction methods analyses, effect sizes are reported as Hedge's g av, calculated using the resources provided by Lakens (2013), and paired t-tests were used to treat each phase-encoding method independently.

Application to G-Ratio Mapping
Our final set of analyses aimed to assess any benefits of the consecutive HySCO WA approach on g-ratio mapping by comparing it to the use of consecutive HySCO AM on the generation of g-ratio maps. As neuroscientific studies utilising the g-ratio are typically performed at the group level, analyses were performed in group space. The g-ratio was calculated from axonal-water fraction (AWF) and MT saturation maps according to Ellerbrock and Mohammadi (2018): with MVF MR being the myelin-volume fraction estimated from the MT saturation map and AVF MR being the axonal-volume fraction. The AVF MR was estimated as AVF MR = (1 − MVF MR ) AWF according to Stikov et al. (2015), where AWF was obtained by combining the intra-cellular fraction (ν icvf ) and isotropic fraction (ν iso ) maps from the neurite orientation dispersion and density imaging (NODDI) toolbox 5 (Zhang et al., 2012) as AWF = (1−ν iso ) ν icvf . The MT saturation map was obtained from the hMRI toolbox as described above. For calibration of the MT saturation map to a myelin-volume fraction map (MVF MR = αMT sat ), we used the g-ratio based calibration method as reported in Ellerbrock and Mohammadi (2018) and Mohammadi and Callaghan (2020), with a reference g-ratio value of 0.7 in the splenium (Mohammadi et al., 2015a), yielding α = 0.217.
Transformation from native to group space was implemented using the voxel-based quantification (VBQ) approach as implemented in the hMRI toolbox (Draganski et al., 2011), which allows for preservation of the quantitative values of the g-ratio maps during smoothing (Tabelow et al., 2019). Briefly, inter-subject registration using DARTEL (Ashburner, 2007) was performed on the segmented MT saturation grey and white matter probability maps. The resulting DARTEL template and deformations were then used to normalize the MT saturation, g-ratio and Jacobian maps to Montreal Neurological Institute (MNI) space at 1.5 mm × 1.5 mm × 1.5 mm resolution with a tissue weighted smoothing kernel of 6mm full width at half maximum.
Differences in g-ratio mapping were determined by comparing the calculated g-ratio values. As the g-ratio can only be defined when the white matter of the diffusion and structural (MT saturation) maps overlap, a g-ratio of zero is recorded for non-overlapping voxels (Mohammadi and Callaghan, 2020). We predicted, therefore, that the poorer overlap between the AM diffusion data and structural white matter map (due to the higher blurring in the AM diffusion data) would increase the number of voxels with a g-ratio of zero per participant, decreasing the averaged g-ratio values at the group-level in regions with high susceptibility-related distortions.
To quantify the reduction of group-averaged g-ratio values within regions suffering from high levels of distortion, for each participant we extracted the mean g-ratio from the AM and WA g-ratio maps from regions previously identified as suffering from distortion. This was performed using a mask created from thresholding the Jacobian map at 20% (see also Comparison of the diffusion MRI phase-encoding combination methods section above). We then compared these to the averaged g-ratio value in white matter using two tailed paired t-tests, thresholded at p < 0.05, calculated in SPSS v25 with effect sizes reported as Hedge's g av, calculated using the resources provided by Lakens (2013). We hypothesised that smaller g-ratio values would be apparent in the AM-based g-ratio map because smoothing with zero-valued g-ratio map voxels, indicative of misalignment, would systematically lower the g-ratio values.
Finally, to localise the voxels where the AM-based g-ratio values were smaller than the WA-based g-ratios, we performed a voxel-wise comparison of the AM and WA g-ratio maps within the same 20% Jacobian distortion mask. We did this using twotailed paired t-tests to identify significant differences in g-ratio values at statistical thresholds of p < 0.05 family wise error (FWE) with a minimum cluster size of 10 voxels.

Comparison of Diffusion MRI Distortion Correction Pipelines
The average root mean square of the difference in FA map alignment between the blip-up and blip-down data was significantly smaller for the consecutive HySCO approach in comparison to the standard HySCO approach for all participants (Figure 3, mean difference = 0.017, t(17) = 8.65, p < 0.001, Hedge's g av = 0.91), indicating better distortion correction by the consecutive HySCO pipeline. Of note, smaller deviations (Figure 3) were observed between the blip-up and blip-down data after FSL distortion correction in comparison to the standard HySCO approach (mean difference = 0.040, t(17) = 10.38, p < 0.001, Hedge's g av = 2.49) and the consecutive HySCO approach (mean difference = 0.023, t(17) = 8.63, p < 0.001, Hedge's g av = 1.72).

Comparison of Diffusion MRI Phase-Encoding Combination Methods
Having established that consecutive HySCO gave better distortion correction than standard HySCO, we then turned to examining the consecutive HySCO AM and consecutive HySCO WA phase encoding combination methods. The improvement to the spatial accuracy gained by using the WA to combine the opposite phase-encoding data is shown in Figure 4. In panels A and B, additional white matter, in line with the anatomically accurate white matter tissue probability derived from the MT saturation map, can be observed in the ventromedial prefrontal cortex when using the WA compared to using the AM. In panels C and D, additional and better definition of the white matter, again in line with the anatomically accurate MT saturation white matter tissue probability map, is evident in the inferior temporal gyrus when using the WA compared to using the AM. Using the WA seems, therefore, to provide better alignment between the diffusion and quantitative structural maps.
To provide a quantitative examination, we also calculated the mean percentage overlap between the diffusion and MT saturation white matter tissue probability maps for the consecutive HySCO AM, consecutive HySCO WA and, as an additional benchmark, FSL (see Figure 5). This quantitative comparison supported the qualitative visual inspection.
Interestingly, across the whole of the white matter map, a larger percentage overlap between the diffusion and MT saturation white matter tissue probability maps was observed when using the AM to combine the FSL distortion corrected blip-up and blip-down data compared to the consecutive HySCO WA (mean difference = 1.64, t(17) = 6.47, p < 0.001, Hedge's g av = 0.88), although there was no difference between the methodologies at 10% distortion (mean difference = 0.34, t(17) = 0.68, p = 0.51, Hedge's g av = 0.13). However, we found that using the WA to combine the consecutive HySCO distortion corrected blip-up and blip-down data resulted in greater percentage overlap of the diffusion and MT saturation white matter tissue probability maps compared to using the AM to combine the FSL distortion corrected data in regions suffering from greater levels of distortion (20% distortion: mean difference = 3.58, t(17) = 6.13, p < 0.001, Hedge's g av = 1.19; 30% distortion: mean difference = 7.26, t(17) = 9.62, p < 0.001, Hedge's g av = 1.93) (Figure 5).
The consecutive HySCO WA method seems, therefore, to be particularly advantageous in regions suffering from high distortion, for example, in the ventromedial prefrontal cortex (as shown in Figures 4A,B), inferior temporal gyri (Figures 4C,D), as well as the middle temporal gyri, medial temporal lobe (including the anterior hippocampi), and the brain stem (Figure 6).

Application to G-Ratio Mapping
Finally, we set out to assess whether there was any benefit of using the consecutive HySCO WA diffusion data over the consecutive HySCO AM data for g-ratio mapping. Figure 7 shows that the coverage of the g-ratio map across white matter is enhanced when using the WA combination compared to the AM combination. The focus is the ventromedial prefrontal cortex and the brain stem, two regions that suffer from high levels of distortion. As can be observed, greater accuracy was obtained in the white matter probability maps and more information was available in the AWF map. Consequently, this resulted in a larger area where the g-ratio could be defined when using the diffusion data corrected using the WA compared to performing the same calculations on the AM combination.
We next sought to examine the g-ratio mapping data at the group level, predicting higher mean g-ratio values when there was better overlap between the diffusion data and the MT saturation white matter tissue probability map. As our previous analyses identified that the WA combination resulted in greater overlap FIGURE 5 | Mean (±standard deviation across participants) percentage overlap of the white matter probability maps created from the different phase-encoding combination diffusion data with the MT saturation white matter tissue probability map. The 10, 20, and 30% distortion graphs focus on regions suffering from greater distortion, determined by thresholding the Jacobian map at distortion values higher than 10, 20, and 30% respectively. HySCO, consecutive HySCO; AM, arithmetic mean; WA, weighted average. **p < 0.01, ***p < 0.001. between the diffusion and MT saturation white matter tissue probability maps in areas suffering from high levels of distortion, we predicted that improvements in g-ratio mapping would be evident when using the WA methodology in such regions.
To investigate this, we first compared the mean g-ratio values from across the whole of white matter, finding significantly greater WA-based g-ratios in comparison to the AM-based g-ratios [WA-based mean = 0.55; AM-based mean = 0.54; mean difference = 0.0071, t(17) = 4.58, p < 0.001, Hedge's g av = 0.45]. To create a non-biased baseline mean g-ratio value, for each participant we calculated a combined mean g-ratio for the AM and WA white matter maps. We then extracted the mean g-ratio from the regions suffering from high levels of distortion (identified by thresholding the Jacobian map at 20%), and subtracted this from the combined mean g-ratio in order to provide a measure of the difference in g-ratio between highly distorted regions and the whole of white matter. We predicted that the difference in mean g-ratio values between the highly distorted regions and the combined mean g-ratio across the whole of the AM and WA white matter maps would be smaller for the WA-based g-ratios. This is because we expected the g-ratio values in the AM-based g-ratio map to be lower in highly distorted regions, and thus further away from the combined mean g-ratio. This would be due to the poorer overlap between the diffusion and MT saturation white matter tissue probability maps increasing the number of voxels with a g-ratio of zero.
Indeed, a smaller difference between the g-ratio values in highly distorted regions and the combined mean g-ratio was apparent for the WA-based g-ratios compared to the AM-based g-ratios (mean difference = 0.092, t(17) = 12.76, p < 0.001, Hedge's g av = 1.98; Figure 8A). This suggests, therefore, that in regions suffering from high levels of distortion, the WA-based g-ratio values were closer to the white-matter average than AM-based g-ratio values, thus providing a better representation of g-ratio values in regions with high susceptibility-related distortions.
Finally, to further localise where, within the regions with high susceptibility-related distortions, the AM-based g-ratios were significantly smaller than the WA-based g-ratios, voxelwise statistics were performed. As can be seen in Figure 8B, higher g-ratio values were identified (top row) when using the WA combination in regions suffering from the greatest levels of FIGURE 6 | Examples of regions where greater distortion was observed across the participants. Images are the average deviation from unity of the Jacobian maps (top row) and average MT saturation map (bottom row) across participants for the whole sample in Montreal Neurological Institute (MNI) space. In the deviation from unity of the Jacobian maps, values away from zero indicate areas suffering from greater distortion (see Eq. 2). Inter-subject registration was performed using DARTEL as implemented in SPM12 (Ashburner, 2007). The DARTEL template and deformations were used to normalize each participant's Jacobians and MT saturation map to the stereotactic space defined by the MNI template (at 1.5 mm × 1.5 mm × 1.5 mm resolution). VMPFC, ventromedial prefrontal cortex.
FIGURE 7 | The impact of reduced blurring artefacts on g-ratio mapping when using the weighted average (WA) combination compared to the arithmetic mean (AM), illustrated for a single participant, in their native space. Depicted in the 1st column are the constituents of the g-ratio: myelin-volume fraction (MVF) based on the MTsat (top), axonal-water fraction (AWF) using AM (middle), and AWF using WA (bottom). In the 2nd column are the associated white matter probability maps (c2): MT c2 (top), AM c2 (middle), and WA c2 (bottom). In the 3rd column is the deviation from unity of the Jacobian map (top, where values away from zero indicate areas suffering from greater distortion, see Eq. 2) and the resulting g-ratio maps AM g-ratio (middle), and WA g-ratio (bottom). The red boxes and crosshairs highlight the ventromedial prefrontal cortex, and the orange boxes the brain stem, both of which suffer from high levels of distortion. Greater accuracy was obtained in the white matter probability maps when using the WA phase-encoding combination method, resulting in a larger area where the g-ratio was defined.
FIGURE 8 | Comparison of g-ratio values when using the weighted average (WA) and arithmetic mean (AM) combinations. (A) Mean (±standard deviation) difference in g-ratios between highly distorted regions (defined by the 20% Jacobian distortion mask) and the whole of white matter. The smaller difference evident after the WA combination suggests the g-ratio values were closer to the white matter average, providing a better representation of g-ratio values in regions with high susceptibility-related distortions. (B) Voxel-wise comparison within the 20% Jacobian distortion mask showing regions with significantly higher g-ratio values when using the WA in comparison to the AM combination (top row), and their overlap with regions suffering from high levels of distortion. This is observable via the deviation from unity of the Jacobian map (bottom row, where values away from zero indicate areas suffering from greater distortion, see Eq. 2). Comparison images are thresholded at p < 0.05 FWE corrected and displayed on the average MT saturation maps for the whole sample. The left images show the ventromedial prefrontal cortex and brain stem, the right images show the inferior and middle temporal gyri. ***p < 0.001. distortion (bottom row), including the ventromedial prefrontal medial frontal cortex, brain stem, and inferior and middle temporal gyri (see also Table 1). No significant voxels were identified for the reverse contrast.

DISCUSSION
Our aims in this study were to improve the susceptibilityrelated distortion correction for diffusion MRI data in SPM, to deploy an alternative approach to reducing image blurring in diffusion MRI data when combining blip-up and blip-down EPI data after susceptibility-related distortion correction, and to assess any benefits of these changes for g-ratio mapping. We found that the new processing pipeline, consecutive HySCO, improved distortion correction. Moreover, the WA compared to the commonly used AM combination method showed a distinct benefit following distortion correction in SPM, and this performed well when benchmarked against FSL, particularly in brain areas with strong susceptibility-related distortions. Finally, we showed that the consecutive HySCO WA performed better than the AM method when combined with other quantitative MRI parameters -MPM -to perform g-ratio mapping. These improvements mean that researchers can conveniently access a wide range of diffusion-related analysis methods within one framework because they are now available within the opensource ACID toolbox as part of SPM, which can be easily combined with other SPM toolboxes, such as the hMRI toolbox, to facilitate computation of myelin biomarkers that are necessary for g-ratio mapping.

The New Diffusion MRI Distortion Correction Pipeline in ACID
Recent work has highlighted that while the distortion correction module in SPM-ACID (HySCO2) performed well in correcting distortion in b0 images, there were limitations in its ability to reduce susceptibility-related distortion (Gu and Eklund, 2019), highlighting a need for improvement. To enhance the performance of the HySCO2 module, we performed two consecutive applications of the HySCO2 module instead of the standard single application, with ECMOCO for motion and eddy current correction performed in both the standard and consecutive HySCO pipelines. The second iteration of HySCO2 in the consecutive HySCO pipeline used FA maps computed after the initial distortion correction as input. We found that the consecutive HySCO pipeline improved distortion correction as compared to the standard HySCO pipeline, while noting that it was less efficient than the FSL TOPUP and eddy pipeline.

WA Combination -A Method to Reduce Susceptibility-Distortion Related Blurring Artefacts
The new consecutive HySCO pipeline was used to test the performance of an alternative phase-encoding combination approach. Areas that are strongly squeezed by susceptibilityrelated distortions suffer from blurring artefacts when taking the AM of the susceptibility-related distortion corrected blip-up and blip-down images. We proposed the WA combination to reduce this blurring artefact. The WA combination uses the Jacobian of the distortion field (estimated by the HySCO2 module in ACID during susceptibility-related distortion correction) to downweight areas that are squeezed in, for example, the blip-up image, while up-weighting the same (stretched) regions in the opposite blip-down image. As illustrated in Figure 4, compared to the consecutive HySCO AM, using the consecutive HySCO WA combination resulted in the identification of additional and better spatially defined white matter, in line with the more anatomically accurate MT saturation white matter tissue probability map that was used here as a reference. The benefits of the WA phase-encoding combination method were also apparent when using quantitative comparisons. Comparing the use of the WA and AM combination methodologies (following distortion, eddy current and motion correction in SPM) identified an increase in the percentage overlap of the diffusion and MT saturation white matter tissue probability maps when using the WA. This was apparent across the whole white matter map, but particularly when focusing on regions suffering from high levels of distortion (e.g., ventromedial prefrontal cortex, inferior and middle temporal gyri, medial temporal lobe, brain stem). The WA also compared well with FSL, which was used for benchmarking, and indeed showed a greater percentage overlap of the diffusion and MT saturation white matter tissue probability maps for regions suffering from high levels of distortion.

Application to G-Ratio Mapping
The beneficial effect of using the consecutive HySCO WA combined diffusion data was illustrated for g-ratio mapping, where blurring effects can reduce the overlap between the diffusion-MRI-based axonal biomarker and the myelin marker maps. The blurring artefacts in diffusion MRI (which are decreased when using a WA combination) lead to a worse delimitation between grey and white matter, which in turn can propagate toward the resulting g-ratio map (Mohammadi and Callaghan, 2020). As the MR g-ratio is defined only in white matter (Stikov et al., 2015), it is typically created by combining an axonal water fraction estimated in the white matter probability maps from diffusion MRI (here, the AWF) and a white matter myelin biomarker map (here, the MT saturation). In current implementations, whenever one of the two constituents is not defined, the g-ratio value cannot be calculated, resulting in undefined areas for the g-ratio map which are set to zero (see Figure 7 in Mohammadi and Callaghan, 2020). By using the consecutive HySCO WA method to combine the diffusion data, higher quality data are available to calculate the AWF, enabling the g-ratio map to cover a larger area of the white matter than when using the consecutive HySCO AM method (illustrated in our Figure 7).
As these undefined areas vary between individuals, we additionally assessed the beneficial effects of the WA over the AM combination for g-ratio mapping at the group level. We predicted that less overlap between the diffusion and MT saturation white matter tissue probability maps would increase the number of voxels with a g-ratio of zero, decreasing the local g-ratio values in the group average. Indeed, this was what we found. In regions suffering from high levels of distortion, the WA-based g-ratio values were closer to the white-matter average than the AM-based g-ratio values. This prediction was further supported by the observation that the AM-based g-ratio values were significantly smaller than WA-based g-ratio values in regions with high susceptibility-related distortions. The improved definition of white matter in regions suffering from high levels of distortion by the WA technique seems, therefore, to also result in better representation of g-ratio values in regions with high susceptibility-related distortions.

Limitations and Future Directions
We developed and tested our consecutive HySCO WA methodology within the SPM framework to facilitate the combination of diffusion and quantitative MT saturation maps, the latter of which has been developed uniquely within SPM (Tabelow et al., 2019). As such, we aimed to improve the tools available in ACID and SPM. It may be that applying the WA methodology developed here to other susceptibility-related distortion techniques, such as FSL, would also provide further benefits, however, examining this was beyond the scope of the current study.
We note that, even with the improvements to the initial SPM distortion correction pipeline, FSL TOPUP and eddy provided better initial distortion, motion and eddy correction than the SPM equivalents. This is potentially due to the consecutive HySCO pipeline using multiple interpolation steps, introducing additional blurring effects. Future work will aim to develop a methodology to combine these interpolation steps, further increasing the effectiveness of HySCO distortion correction. Of note, however, despite the poorer performance of HySCO susceptibility-related distortion correction as compared to FSL, the WA combination in ACID overcame these initial differences in distortion correction between the two toolboxes. Further development of the distortion correction methods in HySCO may result in even greater anatomical accuracy in the future.
The consecutive HySCO WA distortion correction technique may also have implications beyond the demonstrated improvements in g-ratio mapping. For example, it could assist with improved alignment of the diffusion to anatomical data when relating tractography-based macroscopic brain connections to their cortical origin (e.g., Jbabdi et al., 2015), or for multivariate analyses such as those deployed by Draganski et al. (2011). The tools developed here are freely available in the ACID toolbox for SPM for use in other methodological developments such as these.

CONCLUSION
We present a new WA approach to combining blip-up and blip-down EPI images (after susceptibility-related distortion correction) that reduces image blurring in diffusion MRI data. In addition, use of the WA method in combination with other quantitative structural MRI maps that can be processed within SPM improved g-ratio mapping at the group level in regions with high susceptibility-related distortions. This WA approach is now available in the open-source SPM-ACID toolbox together with a new processing pipeline for improved susceptibly-related distortion correction.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by University College London Research Ethics Committee (project ID: 6743/001). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
IC: methodology, investigation, formal analysis, writingoriginal draft, and writing -review and editing. MC: formal analysis, and writing -review and editing. NW: resources and writing -review and editing. EM: conceptualisation, methodology, formal analysis, writing -original draft, writing -review and editing, supervision, and funding acquisition. SM: conceptualisation, methodology, formal analysis, writing -original draft, and writing -review and editing. All authors contributed to the article and approved the submitted version.