Effects of Inaccurate Response Function Calibration on Characteristics of the Fiber Orientation Distribution in Diffusion MRI

Diffusion MRI of the brain enables to quantify white matter fiber orientations noninvasively. Several approaches have been proposed to estimate such characteristics from diffusion MRI data with spherical deconvolution being one of the most widely used methods. Constrained spherical deconvolution requires to define – or derive from the data – a response function, which is used to compute the fiber orientation distribution (FOD). This definition or derivation is not unequivocal and can thus result in different characteristics of the response function which are expected to affect the FOD computation and the subsequent fiber tracking. In this work, we explored the effects of inaccuracies in the shape and scaling factors of the response function on the FOD characteristics. With simulations, we show that underestimation of the shape factor in the response functions has a larger effect on the FOD peaks than overestimation of the shape factor, whereas the latter will cause more spurious peaks. Moreover, crossing fiber populations with a smaller separation angle were more sensitive to the response function inaccuracy than fiber populations with more orthogonal separation angles. Furthermore, the FOD characteristics show deviations as a result of modified shape and scaling factors of the response function. Results with the in vivo data demonstrate that the deviations of the FODs and spurious peaks can further deviate the termination of propagation in fiber tracking. This work highlights the importance of proper definition of the response function and how specific calibration factors can affect the FOD and fiber tractography results.

Diffusion MRI of the brain enables to quantify white matter fiber orientations noninvasively.
3 Several approaches have been proposed to estimate such characteristics from diffusion MRI 4 data with spherical deconvolution being one of the most widely used methods. Constrained 5 spherical deconvolution requires to define -or derive from the data -a response function, 6 which is used to compute the fiber orientation distribution (FOD). This definition or derivation 7 is not unequivocal and can thus result in different characteristics of the response function which 8 are expected to affect the FOD computation and the subsequent fiber tracking. In this work, 9 we explored the effects of inaccuracies in the shape and scaling factors of the response 10 function on the FOD characteristics. With simulations, we show that underestimation of the 11 shape factor in the response functions has a larger effect on the FOD peaks than 12 overestimation of the shape factor, whereas the latter will cause more spurious peaks.

13
Moreover, crossing fiber populations with a smaller separation angle were more sensitive to

26
Diffusion MRI allows to characterize tissue microstructure in vivo and noninvasively by 27 measuring the anisotropic diffusion of water molecules [1,2]. Diffusion tensor imaging (DTI) [3] 28 is the most widely used model in clinical studies to relate the diffusion MRI signals to the 29 diffusion characteristics of the underlying tissue. However, DTI is inadequate to estimate the 30 directional information in voxels containing crossing fibers [4,5]. A commonly used approach 31 to resolve more complex fiber configurations in the brain is spherical deconvolution (SD) [6-32 8]. SD also allows for the extraction of fiber population specific microstructural measures 33 derived from the magnitudes of the fiber orientation distribution (FOD) functions, such as 34 apparent fiber density (AFD) [9] and hindrance modulated orientational anisotropy (HMOA)

36
SD requires an appropriate response function as input to estimate the FOD [7]. The 37 response function, representing the diffusion signal for a single fiber population, is ideally 38 calibrated from the acquired diffusion MRI data [11,12]. In brief, for each subject, the voxels 39 containing only single fiber populations are localized, and an average of the diffusivity 40 characteristics within those voxels is used to represent the subject specific response function.

41
An inadequately chosen response function can affect the quantification of FOD characteristics 42 like AFD and HMOA, as well as the fiber tractography.

43
In order to compare inter-subject AFD, Raffelt and colleagues [9] chose a response 44 function common to all subjects to minimize the differences between subjects for voxel-wise 45 AFD comparison. However, this may potentially result in a bias in the estimated FOD.

46
Specifically, the use of such a common response function for group-wise analysis may cause 47 biases in the FOD peak orientations for individual subjects. Therefore, whereas a common 48 response function is optimal for the comparison of AFD and HMOA in group studies [9], it is 49 unclear whether this is also optimal for group-wise tractography studies because of the 50 potentially inaccurate FOD peak orientations and concomitant spurious FOD peaks. Intuitively, 51 the difference in response function characteristics across healthy subjects are not expected to 52 be large, as response functions are generally averaged from more than hundreds of voxels 53 that are supposed to contain single fiber populations [6,7,12]. This was partly demonstrated by [7,19,24] to improve the conditioning of the deconvolution problem, which is further referred to 127 as constrained SD (i.e., CSD). In addition to directional information, the magnitudes of the FOD The response function used in the CSD process can be either simulated or derived 135 directly from the data. Following the latter approach, which is more common, voxels that have 136 a high chance of containing single fiber populations are used to calibrate the response function.

137
A straightforward approach to numerically implement the concept of a single fiber population 138 is to threshold, for instance, the fractional anisotropy (FA), above a pre-defined value.

139
However, the choice of FA threshold is not trivial and can cause inaccuracies in the response 140 function estimation [12]. A data-driven method using a recursive calibration framework was 141 proposed to estimate the response function from the subject data in an unbiased way [12].

142
This method estimates which voxels contain single fiber populations by iteratively excluding 143 voxels which do not have a single dominant orientation and updating the estimated response 144 function.

145
The choice of the fiber response function has an impact on the peak directions and (1) 167 where and are the axial and the radial diffusivity of the single fiber population, is the ( , )

168
polar angle set between the fiber orientation and the applied gradient. Given the axial 169 symmetry property of the diffusion tensor, Eq. (1) can be simplified as , ( ) = cos 2 + sin 2 = cos 2 + (2) 170 where is the absolute difference between the axial and radial diffusivity. For simplicity, = -

171
if we assume that the signal from each fiber population is a function of , then the ( , ) 172 diffusion-weighted signal can then be rewritten as [3] , ( , ) = 0 where is the non-diffusion-weighted signal and is the b-value that represents the strength = , ( ) = 0 e -( cos 2 + ) 0 e -cos 2 (4) where . Eq. (4) highlights the dependency of on the shape factor and the scaling = e -

177
factor , following the definition in previous studies [8]. In this equation, the scaling factor 178 depends only on the radial diffusivity of the fiber response, representing the isotropic diffusion 179 within the fiber population, whereas the shape factor depends on the difference between the 180 axial and radial diffusivities, representing the anisotropic diffusion within the fiber population. to not alter the shape factor . We can then study the effects of on FOD characteristics, by 192 selectively introducing a discrepancy into the shape or the scale of a simulated single fiber 193 signal with respect to the response function. signal generated by a crossing fiber configuration can then be described by where is the volume fraction of each fiber population, is the total number of fiber 205 populations intercrossing the voxel, and is the angle between the applied gradient and the ( )

206
fiber population. In our work, we focus on configurations of two crossing fiber populations, ℎ 207 but the equations of generating the diffusion-weighted signals can also be extended to analyze 208 the FOD characteristics for more than two fiber populations. were extracted using a Newton-Raphson gradient descent method [28]. All FOD peaks that 217 were smaller than an absolute threshold of 0.1 were regarded as contributions from noise and 218 thus discarded to reduce false positives [29]. All peaks were clustered to the nearest simulated

231
The AFD computation was performed as the integral of the FOD magnitudes assigned 232 to each peak, which in the literature is commonly referred to as "lobe". The calculation of the 233 AFD is similar to what was used in a previous study [30], except that we use the gradients 234 generated by the electromagnetic model [31] to segment the FODs for each fiber population 235 instead of using gradients generated by an icosahedron model.

248
For both simulations, two sets of response functions were tested to achieve (a) different 249 shape but the same scaling factors, by increasing from to 0.6 × 10 -3 mm 2 /s 1. We clustered the peak directions to make sure that we are always comparing the 262 angular deviations between the simulated fiber orientation and the FOD peak orientation most images, and with an isotropic spatial resolution of 1.25 x 1.25 x 1.25 mm 3 . We performed CSD 275 based tractography in ExploreDTI with a step size of 1 mm, an FOD threshold of 0.1, an 276 angular threshold of 30 • , and seeding points per 2mm x 2mm x 2mm across the whole brain.

277
All the tracts were constructed with deterministic fiber tracking to facilitate data interpretation. The reference response function for the in vivo dataset was represented by the diffusion 282 tensor fit to the response function, as estimated with the recursive calibration approach [12].

283
Similar to the method described in Section 2.3.2, the diffusion tensor was used to model the  (1) 301 the region with a single fiber population as identified during the recursive calibration step 302 (referred to as "SFP-mask"); and (2) the region with voxels for which FA > 0.2 (referred to as 303 "FA-mask").

FOD characteristics of single fiber populations
307 308 Fig. 2 shows the effect of changing the shape factor and the scaling factor of the 309 response function on the FOD characteristics in a single fiber population. At SNR < 20, the 310 average number of spurious peaks increases when the shape factor increases, but only slightly 311 increases when the scaling factor decreases ( Fig. 2A). The angular deviation depends mainly 312 on the SNR and is far less affected by changes in shape or scale factor of the response function (Fig. 2B). By contrast, changes in peak magnitude (Fig. 2C) and the AFD (Fig. 2D) as a function 314 of shape and scaling factor of the response function are more pronounced than due to 315 differences in SNR level alone. Notice that the effect of changing the scaling factor (up to 316~60%) is roughly three times larger compared to changing the shape factor (up to ~20%). inter-subject comparisons of AFD studies, which will be discussed further in Section 4.3. At

409
SNR 10, the AFD differences are more related to noise than to the shape of the response 410 function for smaller separation angles (below 60 • ). As for the second peak (Fig. 6B), the AFD 411 can change from -30% to 20% when the shape factor was modified from -50% to 50%, data set. The AFD shows a very different pattern in relation to the shape factor changes 473 compared to scaling factor changes. The AFD differences are homogenous throughout the 474 brain when the scaling factor varies, while the outliers indicate the voxels where there are 475 potential geometrical differences in the estimated AFD from the reference, such as merging or 476 spurious peaks. The AFD differences are up to 98% when the scaling factor decreased by 477 0.2. When changing the shape factor with -0.3 mm 2 /s to 0.3 mm 2 /s, the highest × 10 -3 × 10 -3 478 differences (around 6 to 8%) were observed in areas with a single-fiber population, such as 479 the corpus callosum. Notice that bigger changes of the shape factor α makes the AFD 480 difference more heterogeneous across the brain.  At SNR levels of 30 and 50, the FOD characteristics are consistently affected by the 550 choice of the response functions, while at SNR of 10, noise is the dominating factor that affects 551 the FOD properties (Fig. 3). In addition, more spurious peaks are observed at SNR of 10. At 552 relatively high SNR levels, the shape factor of the response function has a greater impact on 553 the results than the scaling factor. In particular, using a sharper response function for 554 separation angles below 50 • can potentially increase the angular resolution of CSD and can, 555 therefore, improve the estimation of the number of peaks (Fig. 3). The shape of the response 556 function was reported to vary with axonal injury and brain maturation, whereas the scaling 557 factor was observed to change as result of demyelination, axonal diameters and axonal density 558 changes [10,35]. This implies that in brain regions affected by disease, applying CSD with a The extent to which the FODs will be affected by the response function depends largely 565 on the separation angle between crossing fiber populations (Fig. 4). More orthogonally 566 crossing fiber orientations are less sensitive to response function changes, as originally 567 suggested in the spherical deconvolution paper [6]. In voxels containing crossing fiber 568 configurations with smaller separation angles (e.g., below 60 • ), the average angular deviations 569 and their variance increase rapidly with lower shape factors of the response function. By For fiber populations with separation angles below 55 • , CSD fails to estimate the correct 577 number of peaks when response functions with a lower shape factor are employed, leading to 578 artificially higher AFD values (Fig. 6). As FOD peaks merge together when the shape factor is 579 further decreased, the AFD becomes close to the integral of the total FOD amplitudes within 580 the voxel. This is shown in Fig. 6

Effect of FOD angular deviations on fiber tracking
589 If the angular deviations of the FOD peaks are similar in the neighborhood voxels along 591 the white matter pathways, accumulating effects on reconstructed fibers will be significant. By 592 contrast, the heterogonous angular deviations of the FOD peaks may only change the voxel-593 wise characteristics like AFD and number of fiber population peaks, the fiber pathways remains 594 if the angular deviations of FOD was not big enough to end in different voxels in the trajectory.

595
Generally, fiber tractography results will not be severely affected in the main part of the fiber 596 bundles, but may show subtle differences at the edges (Fig. 10). In addition, the termination of 597 fiber pathways passing through crossing regions can be affected [12]. With the in vivo HCP The reference value of the shape and scaling factor of the simulated diffusion-weighted 605 signals match with the values in the corpus callosum as reported before. However, recent 606 studies [37][38][39][40] indicated that the diffusivities of fiber bundles in the brain are not always the 607 same. There is not a full map of diffusivity characteristics of each white matter structure yet.

608
Although our simulation study included the same configurations of crossing fiber bundles in a 609 voxel, in reality, the diffusivities of these crossing fibers may not be identical.

610
In this study, we showed tractography results of an HCP subject using the tensor-fit to