In vivo assessment of anisotropy of apparent magnetic susceptibility in white matter from a single orientation acquisition

Multiple studies have reported a significant dependence of the effective transverse relaxation rate constant (R2*) and the phase of gradient-echo based (GRE) signal on the orientation of white matter fibres in the human brain. It has also been hypothesized that magnetic susceptibility, as obtained by single-orientation quantitative susceptibility mapping (QSM), exhibits such a dependence. In this study, we investigated this hypothesized relationship in a cohort of healthy volunteers. We show that R2* follows the predicted orientation dependence consistently across white matter regions, whereas the apparent magnetic susceptibility is related differently to fibre orientation across the brain and often in a complex non-monotonic manner. In addition, we explored the effect of fractional anisotropy measured by diffusion-weighted MRI on the strength of the orientation dependence and observed only a limited influence in many regions. However, with careful consideration of such an impact and the limitations imposed by the ill-posed nature of the dipole inversion process, it is possible to study magnetic susceptibility anisotropy in specific brain regions with a single orientation acquisition.


Introduction
Previous studies reported an orientation dependence of the gradientecho (GRE) signal in brain white matter (WM), especially for the transverse relaxation rate constant * 2 ( Bender and Klose 2010 ;Lee et al. 2011 ;Rudko et al. 2014 ) and the local resonance frequency offset Δ ( Lee et al. 2010 ;Denk et al. 2011 ). Such an orientation dependence is attributed to the myelin sheath surrounding the axons ( Lee et al. 2012 ;W. Li et al. 2012 ). Wharton and Bowtell (2012) proposed to model a single myelinated fibre as a hollow cylinder (HCFM), to capture the compartmentalisation of proton pools, and to characterise the myelin lipids with a cylindrically symmetric magnetic susceptibility tensor. This model was shown to predict the orientation dependence of both relaxation rate and frequency Abbreviations: ALIC, anterior limb of the internal capsule; AMS, apparent magnetic susceptibility; ATR, anterior thalamic radiation; BCC, body of the corpus callosum; CC, corpus callosum; CG, cingulum; CSF, cerebrospinal fluid; CST, corticospinal tract; FA, fractional anisotropy; GCC, genu of the corpus callosum; GP, globus pallidus, CI, confidence interval; GRE, gradient recalled echo; HCFM, hollow cylinder fibre model; MD, mean diffusivity; MMS, mean magnetic susceptibility; MSA, magnetic susceptibility anisotropy; ODF, orientation distribution function; OR, optic radiation; PLIC, posterior limb of the internal capsule; QSM, quantitative susceptibility mapping; ROI(s), region(s) of interest; SCC, splenium of the corpus callosum; SD, standard deviation. SLF, superior longitudinal fascicle; TA, acquisition time; TE, echo time; TR, repetition time; WM, white matter.
( app ) with fibre orientation can be expressed through the components of the susceptibility tensor of lipid molecules and the fibre-to-field angle , which represents the orientation of the fibre with respect to 0 . Reformulated in terms of the susceptibility tensor of a WM fibre (following the convention of X. , this relationship is as follows where ∥ and ⊥ represent the components of the magnetic susceptibility tensor of the fibre parallel and perpendicular to its orientation, and iso is the sum of ⊥ of the WM fibre and any bulk isotropic susceptibility offset of the tissue (e.g., due to the presence of ironrich ferritin molecules). Appendix A explains the relationship between Equation 1 and the formulation of W.  . Note that in the rest of the paper, ∥ and ⊥ represent the components of the fibre's susceptibility tensor, while the apparent magnetic susceptibility obtained with quantitative susceptibility mapping (QSM, Shmueli et al. 2009 ;de Rochefort et al. 2010 ) is referred to as app .
Ex vivo studies of human corpus callosum ( Lee et al. 2010 ) and porcine optic nerve ( Wharton and Bowtell 2015 ) showed that ⊥ is more diamagnetic than ∥ (i.e., ⊥ < ∥ ), and that the difference ∥ − ⊥ (i.e., magnetic susceptibility anisotropy, MSA) is on the order of 5 to 15 ppb. In vivo studies with multiple head orientations ( W. X. Li et al. 2012 ;Wisnieff et al. 2013 ) reported MSA in the same range of values in most white matter fibre tracts considered. In vivo observations of the orientation dependence of the apparent magnetic susceptibility derived from a single orientation acquisition, on the other hand, are rare. W.  studied AMS in a single ROI in WM in vivo , with the subject turning the head in the scanner, and reported susceptibility anisotropy consistent with the values found via reconstructing the susceptibility tensor.
A different approach was taken by Lancione et al. (2017) , who assessed the orientation dependence of susceptibility across WM in a more practical way using a single orientation acquisition. In their approach, WM susceptibility values at different fibre-to-field angles were pooled from different anatomical brain areas and averaged. It should be noted that, despite the practical benefit, this approach averages variations in MSA across white matter regions, which may obscure the global pattern of orientation dependence. This is especially important since the results of Lancione et al. suggest that ⊥ > ∥ , seemingly at odds with studies that used multiple orientation acquisitions.
Against this background, the purpose of our study was to assess the orientation dependence of AMS in white matter in vivo from a single orientation acquisition while reducing the confounding effect of MSA variation throughout the brain by accounting for the orientation dependence of app within anatomical regions. In the same effort to reduce potential confounding factors, we considered the dependence of the apparent MSA on the fractional anisotropy (FA) of the considered fibres, as Lancione et al. reported that MSA becomes more negative when assessed over fibres with higher FA.

Methods
The cohort of healthy controls included 39 volunteers (age 24-65 yrs., median age 31 yrs., 20 female) without known neurological conditions. The subjects were examined on a 3 T MRI system (Prisma fit, Siemens Healthineers, Erlangen, Germany) using a 20-channel headcoil. The study was approved by the local ethics committee and in accordance with the Declaration of Helsinki. All subjects provided written informed consent.

Data acquisition
To obtain the effective transverse relaxation rate * 2 and the magnetic susceptibility, a 3D proton density weighted (PDw) multiecho GRE sequence was applied using the following parameters: flip angle = 6°; TE 1-5 = 8.12/13.19/19.26/24.33/29.4 ms; TR = 37 Fig. 1. Colour-encoded distributions of fibre-to-field angle in a healthy volunteer in all nine ROIs used in this work as well as the entire WM, excluding crossing fibres (bottom rightmost subfigure). The first echo of the synthetic T1w GRE images was used as anatomical reference. Each subfigure represents the average across a 5 mm thick slab (except for SLF and WM, where the average was taken over a 10 mm thick slab to represent the larger ROIs). The orange arrows represent the direction of 0 . Individual ROIs are shown after erosion as described in the text. ms; monopolar echo readout, bandwidth 1-5 = 280 Hz/px, voxel size = (1 × 1 × 1) mm 3 ; TA = 8 mins 26 s. A second 3D multi-echo GRE sequence with identical parameters but flip angle of 35°served as a T 1 -weighted contrast reference. Two diffusion-weighted echo planar imaging scans (voxel size = (1.5 × 1.5 × 1.5) mm 3 , TR = 3300 ms, TE = 84 ms, TA = 6 mins each) with multi-band readout, reversed phaseencoding polarities, and multi-shell diffusion scheme (8, 16, 32, and 48 directions with b-values of 0, 835, 1665 and 2500 s/mm 2 , respectively) were collected to determine the orientation of WM fibre bundles with respect to the static magnetic field, 0 .

Diffusion MRI
Diffusion weighted acquisitions were combined and distortioncorrected using FSL (topup version 5.0.9, eddy_openmp version 5.0.8; Andersson et al. 2003 ;Andersson and Sotiropoulos 2016 ). The diffusion tensor model was fitted to obtain fractional anisotropy (FA) and mean diffusivity (MD) using Mrtrix3 (version 3.0; Tournier et al. 2019 ). The predominant fibre orientation was determined as the angle between the first peak of the orientation distribution function (ODF) and the 0 direction using constrained spherical deconvolution ( Tournier et al. 2007 ). MD maps were rigidly registered to the first echo images of the T1w GRE sequence using Aladin from NiftyReg ( Ourselin et al. 2001 ), as distributed with nicMSlesions (release candidate 0.2; Valverde et al. 2018 ). FA and maps were transformed to the multi-echo T1w space based on this registration.

Segmentation and region-of-interest definition
Using Freesurfer (version 6.0.0; Desikan et al. 2006 ), white matter segmentation was obtained from a synthetic T1w contrast (generated from the two GRE acquisitions with mri_synthesize, TR = 20 ms, FA = 30°, TE = 5 ms), and five-tissue-type segmentation from Mrtrix3 ( Smith et al. 2012 ). The resulting WM mask was additionally eroded with a sphere of radius 2 mm to avoid the effect of potential segmentation errors. To exclude crossing fibres, voxels where the second peak of their ODF was greater than 30% of the first peak were removed from the WM mask (see Figure 1 , right bottom image). All further analyses were restricted to the data from this mask. Additional selections based on FA values are discussed below. The white matter mask described above was segmented into ROIs based on the diffusion data with Tractseg (version 1.7.1; Wasserthal et al. 2018 ). Resulting tract probability maps were transformed to the GRE space, limited to probability > 90% and additionally eroded with a 6-connected structural element. ROIs were selected on the basis of relatively homogeneous fibre orientation, sufficient proportion of voxels with high FA, and in an attempt to minimise their spatial overlap (see Figure 1 , corresponding images). In particular, the following ROIs were considered: corpus callosum (CC), subdivided into the genu (GCC), the body (BCC), and the splenium (SCC); anterior thalamic radiation (ATR), subdivided into the anterior limb of internal capsule (ALIC) and the remaining prefrontal part (further referred to as ATR); corticospinal tract (CST); superior longitudinal fascicles I, II, and III combined in one single ROI indicated by SLF; cingulum (CG); and optic radiation (OR). Since the definition of the corpus callosum in Tractseg included the entire length of commissural fibres crossing CC, we defined CC as the intersection between Tractseg based fibres and dilated Freesurfer labels (see appendix B for more details). Furthermore, the anterior limb of internal capsule (ALIC) was separated from the rest of ATR using the mask of subcortical grey matter obtained from Freesurfer segmentation, and additionally dilated with a spherical element of radius 8 mm. While the posterior limb of the internal capsule (PLIC) can be similarly separated from the corticospinal tract, the typical range of and FA found in PLIC was too narrow for the analysis to justify such subdivision.

Gradient-echo data processing
For each subject, the PDw images of the low flip angle GRE sequence were rigidly registered to the T1w GRE data using Aladin. * 2 maps were calculated from the PDw GRE data via voxel-wise monoexponential linear fitting. Phase images of the individual echoes of the multi-echo GRE data were independently unwrapped using a 3D pathfollowing algorithm ( Herráez et al. 2002 ) from scikit-image (version 0.16.1; Van der Walt et al. 2014 ). Since spatially unwrapped phase of individual echoes can still be affected by arbitrary global multiples of 2 -offsets between individual echoes, these offsets were estimated from the temporal evolution of the mean phase within a sufficiently large and homogeneous ROI (the WM mask) ( Schweser et al. 2016 ). To estimate the time-independent phase offset, a voxel-wise linear fit was performed. The resulting phase offset map was further smoothed by fitting a 3D polynomial of 5th degree and subsequently subtracted from the unwrapped phase of all echoes ( Schweser et al. 2011 ). This corrected unwrapped phase was divided by ⋅ 0 ⋅ T E ⋅ 10 −9 for each echo, where T E is the corresponding echo time, to produce maps of relative difference field in ppb ( Kee et al. 2017 ). These maps were averaged over the last three echoes (TE of 19-29 ms). A phase reliability mask was produced that met the following three criteria: the local SNR estimate of the magnitude image was greater than 15, all components of the phase gradients in all three directions were less than ∕2 by absolute value, and standard deviations of the relative difference field estimated from the three last echoes were less than 4 ppb. The mask was morphologically closed and filled for subsequent inclusion in the background field removal with V-SHARP ( Schweser et al. 2011 ;Wu et al. 2012 ) (spherical kernel with radii 3-8 mm; in-house implementation). QSM values were obtained with iLSQR (STI Suite, version 3.0; Li et al. 2015 ) using the reliability mask described above. QSM values were referenced to the average magnetic susceptibility of CSF in the lateral ventricles (the choroid plexus was removed by segmentation based on the * 2 map according to the procedure outlined in appendix C ).

Data analysis
For further analysis, values of FA, , * 2 , and QSM were extracted from the aforementioned individual tract ROIs for each subject. Unless explicitly noted that the analysis was performed on an individual subject basis, all voxels were considered for analysis regardless of which volunteer they came from (group-averaged approach).
Two alternative approaches were adopted for the additional selection of voxels used to assess the orientation dependence of * 2 and AMS. (1) All voxels with FA > 0.6 were used to assess both * 2 ( ) and AMS( ) without further consideration of FA, or (2) voxels with FA > 0.3 were grouped in four equally populated bins based on the quartiles of the FA distribution in an ROI and independently processed for each bin (for AMS only). In both approaches, selected voxels were further grouped into 10 equally populated bins based on the quantiles of the fibre-tofield angle distribution ( ( 0 . 0 , 0 . 1 ]; ( 0 . 1 , 0 . 2 ]; … ; ( 0 . 9 , 1 . 0 ] ) for visualization purposes. The means and standard deviations were calculated in each bin for AMS (and * 2 ). To assess the orientation dependence of AMS quantitatively, Equation 1 was fitted to the corresponding 10 mean QSM values. Regression was performed with linear least squares (statsmodels 0.12.0, Seabold and Perktold 2010 ) using co s 2 of the mean value in each bin as the independent variable. Note that such an approach averages across multiple voxels pooled from different spatial locations within a given tract ROI. To highlight the difference between such an approach and the voxel-wise analysis of a single specimen rotated in the scanner, we introduce instead of ∥ − ⊥ in Equation 1 and refer to it as the apparent magnetic susceptibility anisotropy, since the obtained estimate may not reflect the difference between the components of any susceptibility tensor. For completeness, iso is referred to as the isotropic apparent magnetic susceptibility. Figure 2 shows the group-averaged distribution of fibre-to-field angles in each ROI as well as for all voxels taken together (see subfigure on the bottom right). Since the ROIs are based on anatomical WM tracts, most of them span only a limited range of the fibre orientation angle, which is also reflected in the -histograms. Taken together, the ROIs are a good representation for the entire range of orientations, since the distribution of angles is skewed towards 90°, as expected, due to purely geometrical considerations (see last subfigure in Figure 2 ). In particular, in the hypothetical scenario where the fibre orientations are uniformly distributed over the solid angle 4 , the number of fibres oriented in the angle interval to + is proportional to sin ( Denk et al. 2011 ). Some of the ROIs selected for this study, however, showed prominent disparity in the number of voxels available for analysis. Specifically, SLF combines three relatively large tracts, whereas ALIC and the rest of ATR comprise an order of magnitude fewer voxels (on average 700 voxels per subject for ALIC and prefrontal ATR compared to 9000 voxels for SLF). Such differences in the number of voxels are reflected in the confidence intervals associated with the estimated model parameters.

Orientation dependence of * 2 in multiple ROIs
Figure 3 displays the orientation dependence of the effective transverse relaxation rate constant * 2 for the different WM ROIs. For most ROIs, * 2 exhibits a very consistent pattern of increasing values with , which matches with observations from previous studies ( Bender and Klose 2010 ;Denk et al. 2011 ;Kor et al. 2019 ;Lee et al. 2017 ;Lee et al. 2011 ). A deviation from this expected qualitative behaviour is observed only in the body of corpus callosum (BCC), where the average * 2 appears to decrease for > 75 • -an angle interval in which the variation of * 2 increases for the entire corpus callosum (CC). Moreover, as can be seen from the overlay subfigure in Figure 3 , * 2 ( ) agrees quantitatively in all ROIs except for BCC (blue dotted line). On average * 2 in BCC is decreased by 1.2 s − 1 , i.e., by about 40% of the full Δ * 2 span.

Orientation dependence of AMS in multiple ROIs
Figure 4 displays the dependence of the apparent magnetic susceptibility app in each considered WM ROI along with the fitted model ( Equation 1 ). In contrast to * 2 , the apparent magnetic susceptibility in these ROIs shows quite different patterns of orientation dependence, which can be summarized as follows: • The splenium of the corpus callosum (SCC) and the optic radiation (OR) follow approximately the expected co s 2 pattern (coefficient of    exhibit a strong non-monotonic orientation dependence of mean AMS with the minimum around = 70 • , which is reflected in the relatively low fit quality ( 2 = 0 . 18 and 0 . 53 , = 0 . 2 and 0 . 02 for GCC and BCC, respectively). The 95% CI for the genu covers a wide range of [-53, 14] ppb, while CI for BCC is strictly positive and spans [8, 62] ppb.
A summary of the estimated regression coefficients along with the 95% confidence intervals and coefficients of determination 2 can be found in Figure 5 and Supplementary Table S1 (including p -values).

Estimation of susceptibility parameters on an individual subject basis
The results presented so far were obtained from a group-averaged analysis and captured the patterns with the highest SNR achievable in this study. The same analysis was repeated at the subject level to assess the amount of variance that can be expected from individual acquisitions. Supplementary Figure S1 shows the distribution of the apparent magnetic susceptibility anisotropy and the isotropic apparent magnetic susceptibility iso (see Equation 1 ) for the individual subjects together with the group-averaged estimates.
For each of the five ROIs with significant orientation dependence in the above group analysis, more than half of the subjects showed significant orientation dependence (at the = 5% level), as indicated by the distributions of the -values of the individual regressions presented in Supplementary Figure S2.

The effect of fractional anisotropy on the orientation dependence of AMS
To assess the role of fractional anisotropy, the analysis presented above was repeated for voxels from different FA groups. In particular, voxels in each ROI were grouped into four equally populated bins using FA quartiles. Figure 6 illustrates this grouping on the joint distributions of fibre-to-field angle and fractional anisotropy (includes only values of FA > 0.3). Equation 1 was then fitted to apparent magnetic susceptibility within the four FA bins (following the second approach described in the Section 2.3 ), yielding four iso and estimates for each WM ROI. Of note to Figure 6 , the joint distributions in GCC, BCC, OR, and to some extent CG exhibit weak correlation. This implies that the different bins used in the analysis above contain voxels that have different FA distributions (further discussed in Section 4.5 .). Figure 7 shows the results of the apparent MSA for the different FA quartile bins. In the following, we summarize the salient observations for the different ROIs.
• Optic radiation (OR). As FA increases, the orientation dependence becomes more prominent, as reflected in the bar plot showing increasing apparent MSA. Note that the values estimated from the first, second, and the last FA bin do not fall within the confidence interval of the reference estimate, indicating a substantial monotonic change in with FA. The coefficient of determination of each regression also increases with increasing FA, indicating better goodness of fit (see Supplementary Figure S3).
• Splenium of the corpus callosum (SCC). All estimates fall within the reference confidence interval, indicating no substantial correlation between and FA. At the same time, the available range of decreases with increasing FA, which is reflected in the increasing width of the CI for the apparent MSA estimates and the concomitant decrease of 2 . The overall range of FA values observed in SCC is also narrower than in OR (see also Figure 6 ).
The apparent MSA increases slightly (i.e., within the confidence interval of the reference result) with increasing FA. 2 is relatively low and 0.5 at most. • Cingulum (CG). With the reference apparent MSA being virtually zero, there is a slight increase in with FA; yet the orientation dependence is inconclusive and 2 is low in all FA bins.
• Superior longitudinal fascicles ( SLF). The apparent MSA decreases (increases by the absolute value) within the reference CI, while 2 stays high for all FA bins. • Corticospinal tract (CST). Both apparent MSA and 2 increase by the absolute value with increasing FA. • Genu of corpus callosum (GCC). Similar to OR, GCC exhibits a correlation between and FA (see Figure 6 ), with most of high FA voxels spanning a narrow range of . With increasing FA, the apparent MSA increases by absolute value, extending outside of the reference CI. • Body of corpus callosum (BCC). This ROI exhibits a complex joint ( , FA) distribution (see Figure 6 ). Shifting the FA window affects the selection substantially. The changes in range from [25, 60] ppb to [-65, 15] ppb.

Discussion
In this study, we examined the orientation dependence of both apparent * 2 and apparent magnetic susceptibility in individual tract-based regions-of-interest in a cohort of 39 healthy volunteers. To counteract the reduction in the number of available voxels caused by the subdivision of white matter into individual ROIs, we first assessed the pattern of orientation dependence in the data averaged across all subjects. Subsequent comparison of such group-average and per-subject analyses revealed that the former reliably represented the median behaviour across the cohort of volunteers. Dashed lines represent quartiles of (marginal) FA distributions within each ROI. Solid black line indicates FA = 0.6 (i.e., the threshold used to select the data in all of the analyses presented above).  Table S1). The light-grey shaded region around the baseline indicates the 95% CI (see Supplementary Table S1 and Figure 5 ). The error bars represent the 95% CI of the corresponding apparent MSA estimates in each FA bin. The variable widths of the dark-grey bars reflect the different widths of the FA bins.

Orientation dependence patterns of * 2 and AMS
The apparent transverse relaxation rate constant * 2 shows a pattern of orientation dependence that is consistent for most of the WM ROIs and is also consistent with the results of previous studies ( Denk et al. 2011 ;Chen et al. 2013 ;Wharton and Bowtell 2013 ;Kor et al. 2019 ). In contrast, the apparent magnetic susceptibility is characterized by distinctly more heterogeneous patterns of the orientation dependence as well as greater scatter between subjects. The latter can be partially explained by the referencing of the QSM values to the susceptibility of CSF, which introduces additional uncertainty in the offset of the individual curves (see Figure 4 ) due to possible segmentation errors of the lateral ventricles as well as the statistical error of QSM values.
In accordance with previously reported observations ( Lee et al. 2010 ;W. Li et al. 2012 ;Wharton and Bowtell 2012 ;Wharton and Bowtell 2015 ), magnetic susceptibility in OR and SCC becomes more diamagnetic with increasing fibre-to-field angle (positive ). In contrast, it becomes less diamagnetic in SLF and CST with increasing fibre angle (negative ). Finally, both genu and body of the corpus callosum (BCC) as well as the anterior limb of internal capsule (ALIC) exhibit a nonmonotonic orientation dependence.
From a biophysical perspective, the observed positive susceptibility anisotropy of WM fibres is thought to arise primarily from the highly structured arrangement of fatty-acid chains in myelin phospholipids ( W. van Gelderen et al. 2015 ), which are themselves characterised by anisotropic magnetic susceptibility ( Lounila et al. 1994 ). The effect of other myelin components such as cholesterol or membrane proteins has not been fully elucidated, although to a limited extent they are thought not to interfere with the resulting co s 2 dependence ( W. . The possibility of explaining the presumed sign reversal of the apparent MSA by the complex myelin content and its variation across WM ROIs remains an open question that is beyond the scope of this work. We go on to discuss each of the empirically observed patterns of AMS orientation dependence and place them in the context of previously reported experimental results.

Orientation dependence characterized by positive MSA
The orientation dependence of AMS in optic radiation and SCC follows the empirical model of Equation 1 quite well over the available range of ( 2 = 0 . 86 in both ROIs). It also agrees with previously reported observations of WM fibres that appear most diamagnetic when oriented perpendicular to 0 . This allows the numerical values of the apparent magnetic susceptibility anisotropy obtained in this study to be related to the magnetic susceptibility anisotropy of myelin lipids as well as their volume fraction (see Equation A1 ), potentially allowing inferences to be drawn about the state of myelin from an assessment of the apparent MSA, , in these ROIs.

Orientation dependence characterized by negative MSA
The pattern of increasing apparent magnetic susceptibility with the fibre-to-field angle, as observed in SLF and CST, contradicts the expectations just discussed. Remarkably, however, Lancione et al. (2017) reported negative MSA values using a similar approach of averaging QSM values across voxels with contributions from potentially very different locations in WM. In addition, X.  reported that magnetic susceptibility anisotropy retrieved from an acquisition with multiple head orientations appeared negative in the posterior limb of internal capsule (PLIC), which constitutes the most anisotropic part of the corti-cospinal tract (according to the definition used in this study). Based on the prevalence of positive MSA in other white matter regions, Li et al. concluded that the negative susceptibility anisotropy observed in PLIC might be artificial, resulting from incomplete background field removal or erroneous field inversion due to the close proximity of the strongly paramagnetic globus pallidus (GP).
Notably, the susceptibility maps obtained in our study seem to be affected by a common inhomogeneity near the GP, as seen in Supplementary Figure S4. Such an artefact appears to decrease the susceptibility immediately posterior to the GP (despite the predominantly parallel orientation of the fibres in the PLIC), while superior to the GP it manifests as a patch of increased AMS vaguely propagating in the direction of 0 (despite the predominantly perpendicular orientation of the ALIC). Importantly, such an overlap of a fibre tract with an artefact caused by another adjacent anatomical structure may produce spurious correlations in the apparent orientation dependence of QSM values. This spatial arrangement, together with the non-monotonic dependence of AMS in ALIC, suggests that the estimate of apparent MSA as well as the isotropic AMS in both limbs of the internal capsule are likely biased and should be viewed with caution. Note that such co-localization of WM fibres with dipole inversion artefacts poses a challenge to any approach that attempts to reconstruct an orientation dependence of AMS by populating different orientations with voxels from different spatial locations. Partitioning WM into individual ROIs reduces the possibility of averaging over such localised perturbation contributions, resulting in less biased estimates of MSA in ROIs without strong artefacts (such as OR).
Although the data in the CST appeared to be contaminated by the artefacts in QSM, the estimated values of apparent MSA were very close to those obtained in SLF, which showed no evident artificial susceptibility variations. It should also be noted that the orientation dependence of AMS in SLF was the most robust in our cohort, with 38 of 39 subjects exhibiting negative and significant apparent MSA. This robustness and the size of SLF suggest that it may be one of the most important contributors to the overall negative apparent MSA found by Lancione et al. (2017) . In addition, the fact that the ROI with the most extensive fibreto-field angle coverage used in the present study showed negative apparent MSA may point to an important limitation of the core assumption of the work. Namely, that aggregation of the reconstructed QSM over multiple voxels with different fibre orientations can essentially substitute for susceptibility tensor reconstruction from multiple orientations in a given voxel. The hypothesis of the artificial origin of negative apparent MSA finds additional support in comparison with the results of Wisnieff et al. (2013) , discussed later in Section 4.3 .

Non-monotonic orientation dependence in corpus callosum
The non-monotonic behaviour of the apparent magnetic susceptibility in the BCC may have its origin in a rather pronounced pattern of darkening on the QSM maps, which is particularly well seen in sagittal view (see the white arrow in Supplementary Figure S5) and appears to propagate into the BCC from the border to the lateral ventricles and CSF. Most notably, the darkening does not appear to correlate with the changes in the fibre-to-field angle and exhibits sharp boundaries in the regions of homogeneous fibre orientation. Moreover, in different subjects, this inhomogeneity tends to occur in different parts of the BCC, possibly in response to head tilt (note that such tilt generally does not affect the distribution in BCC, as fibres remain perpendicular to 0 ). In some subjects, a shift in inhomogeneity was observed towards an area where the boundary between CC and lateral ventricles is oriented mainly perpendicular to the B 0 direction (compare the QSMs on the left and the right in Supplementary Figure S5).
These observations suggest the artificial nature of the pattern. One explanation could be the inability of the scalar dipole model used in QSM to faithfully represent the field produced by strongly anisotropic fibres oriented perpendicular to B 0 and juxtaposed to CSF, which is char-acterised by a different magnetic susceptibility ( Wharton and Bowtell 2015 ;Liu 2019 ).

On the relationship between the orientation dependence of * 2 and AMS
Having discussed the orientation dependence of AMS in individual ROIs, we now turn to the relationship between this orientation dependence and * 2 ( ) . Lee et al. (2011) proposed that the orientation dependence of * 2 in parallel cylinders oriented at an angle to 0 can be related to their magnetic susceptibility as * 2 ( ) = 0 + 1 sin 2 , where itself can be orientation dependent, equivalent to app in Equation 1 , while 0 and 1 are assumed to be constant. Importantly, Lee and colleagues' work did not consider the relationship between the microscopic tensor susceptibility and the observed apparent susceptibility, but instead assumed that * 2 ( ) is driven by ( ) ∼ sin 2 , which allowed them to derive an empirical relationship between the effective relaxation rate and the fibre angle .
Excluding the ROIs for which the orientation dependence of AMS is considered to be affected by one of the artefacts discussed above, we are left with at least OR, SCC, and SLF, which show significantly different patterns of orientation dependence for AMS but not for * 2 . Indeed, the apparent * 2 ( ) in SLF agrees very well with that in SCC, while the data from OR are well within the error margins of the other two ROIs. Such agreement of * 2 ( ) and disparity of app ( ) does not seem consistent with the model of Equation 2 . We attribute such discrepancy to the interpretation of from Equation 2 as the orientation dependent app . A more thorough analysis was conducted by Wharton and Bowtell (2013) , who proposed that * 2 follows a sin 4 dependence, weighted by a quadratic polynomial of MSA and the mean magnetic susceptibility (MMS, see Equation A4 ). In principle, such a dependence allows the sign of MSA to change without substantial qualitative changes in * 2 ( ) , although the quantitative estimation of such an effect exceeds the scope of this work.

Assessment of the obtained numerical values of apparent MSA
Previously reported values for magnetic susceptibility anisotropy of white matter tend to agree that MSA is positive (i.e., ∥ − ⊥ > 0 ) and that white matter fibres appear more diamagnetic when perpendicular to the direction of the magnetic field 0 . This has been demonstrated, for example, with measurements in human corpus callosum ex vivo at 7 T ( Lee et al. 2010 ), human WM in vivo at 7 T ( Wharton and Bowtell 2012 ), porcine optic nerve at 7 T ( Wharton and Bowtell 2015 ), human occipital WM in vivo at 3 T ( W. , and in human spinal cord ex vivo with MR-independent torque measurements at 7 T ( van Gelderen et al. 2015 ). The results indicate the positive sign as well as the range of values of the magnetic susceptibility anisotropy in WM fibres. In our study, the WM ROIs BCC, OR, SCC, and ALIC all show positive apparent MSA, in agreement with the literature. However, the observed susceptibility anisotropy of 20 to 35 ppb appears to be larger than previously reported values (5-20 ppb). This discrepancy is likely due to differences in the QSM processing pipelines used. In particular, the choice of parameters used in the background field removal step (in this work, the radii of V-SHARP convolution kernels and the strength of applied regularisation) affects the overall range of susceptibility values obtained from the dipole inversion. Consequently, the scale of the numerical values of the estimated parameters is also affected, but crucially, the effect is linear (with respect to the V-SHARP regularisation parameter) and identical across all ROIs.
In vivo MSA evaluation using multiple head orientations was reported in volunteers by X.  at 7 T and by Wisnieff et al. (2013) at 3 T. Both studies derived MSA from a cylindrically symmetrical susceptibility tensor. Figure 8 compares the results between our current study Fig. 8. Comparison between the mean apparent MSA ( ) reported in the current study (blue) and the MSA ( ∥ − ⊥ ) found in the other two independent multiorientation studies. The error bars for the present study represent the 95% CI for the group-averaged estimate; in the study of Wisnieff et al. they indicate the standard deviation (SD) in individual subjects; whereas in the study by X. Li et al. they represent the SD of the mean estimates across five subjects. Note that in the latter study the values were reported for posterior thalamic radiation, which overlaps substantially with optic radiation and is therefore labelled OR here. and these two multi-orientation studies in the ROIs common to all three studies.
First, it should be noted that the values of apparent MSA in our study are comparable to those estimated by Wisnieff et al. from acquisitions with three different head orientations, especially for the optic radiation (OR) and the body of corpus callosum (BCC). Furthermore, the estimates of MSA in the splenium (SCC) reported in that study show greater variability compared with those in OR, with one subject having a negative mean MSA. This is consistent with the distribution of apparent MSA estimated from individual subjects in our study (see Supplementary Figure  S1). In the genu of corpus callosum (GCC), the mean MSA estimated in the three-orientation study by Wisnieff et al. was found negative for all four subjects, highlighting the challenge of reliably estimating susceptibility anisotropy in this ROI. Interestingly, using a higher number of head orientations yielded positive MSA in GCC. The authors attributed such a pattern to the difficulty of field inversion (in both the scalar and tensor approaches) near the frontal sinuses. We believe that the same reasoning can explain the negative apparent MSA obtained in GCC in our study and indirectly supports the idea of an artificial nature of the negative in SLF.

Distribution of values on a subject basis
Partitioning WM into individual ROIs and further filtering out voxels based on their fractional anisotropy leads to substantial challenges in estimating the apparent MSA in individual subjects. The distributions presented in Supplementary Figure S1 give an idea of the variations to be expected from datasets acquired at 3 T within clinically acceptable scan times. Only the optic radiation, superior longitudinal fascicles, and the corticospinal tract show an orientation dependence robust enough to yield the same sign of apparent MSA for all subjects. However, even for these ROIs, the average individual 95% confidence interval covers 50% to 70% of the full range of apparent MSA observed in our cohort. Such high uncertainty of individual estimates limits the ability to characterise potential underlying microstructural changes on an individual basis.

Role of fractional anisotropy
Following the work by Lancione et al. (2017) , who showed that the AMS pooled from the entire brain exhibits stronger orientation dependence in voxels with higher FA, we divided the data from each ROI into four equally populated FA bins and compared the estimated apparent MSA. In 5 out of 9 ROIs (namely, SCC, CG, SLF, ATR, and ALIC), δχ estimated from the different FA ranges did not change significantly compared with the reference MSA obtained from the broad window of FA > 0.6. In contrast, the optic radiation showed a monotonic increase in apparent MSA and R 2 with increasing FA. In particular, δχ changed from 9 ppb with FA ∈ ( 0 . 30 , 0 . 47 ] to 46 ppb with FA ∈ ( 0 . 77 , 0 . 86 ] . In turn, for GCC and BCC, changing the FA ranges of the fibres resulted in substantial changes in orientation dependence of AMS. This behaviour may indicate a systematic error, such as the confounding artificial susceptibility inhomogeneities discussed earlier. We are not aware of any model that proposes to explicitly use FA to explain the orientation dependence of QSM values in WM; nevertheless, it is conceivable that the anisotropy of subvoxel distributions of white matter fibres may influence the apparent orientation dependence of quantities derived from measured GRE signals. For example, Kaden et al. (2020) formulated an expression for the macroscopically observed signal shift based on an intravoxel orientation heterogeneity. Similarly, Wharton and Bowtell (2013) proposed a model for the orientation dependence of * 2 and suggested to account for the dispersion of fibre directions within a voxel using a factor λ, analogous to FA. Note that the model proposed by W.  explicitly refers to the myelin lipid volume fraction (see Equation A1 ; also Schweser et al. 2012 ). Lee et al. (2017) proposed an empirical model for * 2 orientation dependence in white matter that successfully used myelin water fraction (MWF) as a proxy for myelin concentration (and thus myelin volume fraction). Uddin et al. (2019) reported a positive correlation between FA and MWF across multiple brain structures when taken together, with the strength of the correlation varying between ROIs. All of this suggests that interpretation of changes in apparent MSA requires a model of AMS orientation dependence that includes both a robust measure of myelin concentration and an accurate model of the influence of FA on apparent magnetic susceptibility. In the absence of such a model, caution should be exercised when interpreting AMS orientation dependence in regions with correlated FA and distributions (see BCC and OR in Figure 6 ).

Limitations to the observed orientation dependence of the AMS
While AMS in OR and SCC seems to follow the phenomenological model for app (θ) ( Equation 1 ), the observations in CC, SLF, CG, and ATR indicate limitations of either the model or the data collection and processing. Therefore, it is important to list these explicitly.
First, apparent magnetic susceptibility as recovered by QSM assumes scalar magnetic susceptibility and is inconsistent with the model of signal formation from a fibre bundle with anisotropic susceptibility. Wharton and Bowtell showed that the perturbation field ΔB z caused by anisotropic susceptibility includes contributions from x -and y-components of the magnetization vector ( Wharton and Bowtell 2012 ). These contributions are not consistent with the scalar dipole deconvolution model, and depend on the orientation of the fibre. Consequently, they can cause AMS to deviate in complex ways from co s 2 θ.
Second, the hollow cylinder model implies a non-linear evolution of the phase of the GRE signal from a white matter fibre ( Wharton and Bowtell 2012 ). Such nonlinearities cause the frequency estimates to depend on echo-time (TE), which has been shown to affect QSM values ( Sood et al. 2017 ;Cronin et al. 2017 ). While the hollow cylinder model predicts both temporal and orientation dependence of the average frequency shift in WM, the model for the orientation dependence in apparent magnetic susceptibility ( Equation 1 ) does not address any TE dependence of QSM values. The data collected in this study predominantly exhibit TE dependence in the first three echoes (8-19 ms), which were intentionally omitted from this work.
Third, more sophisticated theoretical models of the origin of the frequency shift in the GRE signal, such as the Generalized Lorentzian Tensor Approach, GLTA, ( Yablonskiy and Sukstanskii 2018 ), suggest more complex ways to relate the frequency of the observed signal to the underlying model of WM microstructure, implying other sources of inaccuracy in the simplified empirical model of AMS orientation dependence.
All of these factors may therefore help explain the inconsistent AMS orientation dependence observed across different ROIs in this study.

Conclusion
Previously reported differences in magnetic susceptibility of WM fibres oriented perpendicular or parallel to 0 were investigated in a large cohort of volunteers in vivo in a series of WM tracts based on single orientation acquisitions. The expected decrease in AMS with increasing fibreto-field angle was observed in both the splenium of the corpus callosum and optic radiation and tended to agree qualitatively and to some extent also quantitatively with previous observations. Since such orientation dependence is thought to result from the properties and arrangement of myelin sheaths, its assessment could open a way to study microstructural changes in white matter, for example, in the course of a neurodegenerative disease such as multiple sclerosis. However, AMS behaviour in other WM brain regions of interest was found to range from insignificant (cingulum and prefrontal anterior thalamic radiation) to contradictory to studies using multiple head orientations (superior longitudinal fascicles, the body and genu of corpus callosum, the anterior limb of internal capsule). This suggests that the tensor nature of magnetic susceptibility in white matter fibres and multi-exponential evolution of the GRE signal substantially limit the ability of apparent magnetic susceptibility to reflect the macroscopic structure of WM in general. Moreover, presented in this work analysis of the orientation dependence of AMS in individual ROIs appears sensitive to systematic errors and artefacts caused by the ill-posed nature of QSM and the possible inconsistencies between the scalar susceptibility model and the observed phase data. Since the orientation dependence of * 2 in WM has been previously associated with the anisotropy of the magnetic susceptibility tensor, the reported discrepancy between the observed orientation dependences of AMS and * 2 ( ) supports the conclusion of a complex indirect relationship between the apparent magnetic susceptibility and the underlying microstructure. Finally, the fractional anisotropy of WM fibres shows a limited and not yet fully understood correlation with the apparent strength of the orientation dependence of AMS in optic radiation. Based on the above, we expect the most interesting results from further investigations of the possibility to infer microstructure parameters (such as the components of susceptibility tensor) from the evolution of the GRE signal using more sophisticated models.
Despite the above limitations, the strong agreement of the orientation dependence of QSM values observed in the optic radiation and the splenium with the empirical model suggests that a simplified model can be used to extract physically meaningful information in a limited selection of white matter regions. However, assessing the limits of applicability of such an empirical model requires further development in realistic models of GRE signal formation that consistently describe and incorporate other microstructural parameters, such as, for example, fractional anisotropy.

Data and code availability statement
All data used in this study were newly acquired by the authors and present private clinical data of the investigated individuals, and as such cannot be shared or processed outside of the institution where they were acquired due to institutional regulation on personal information privacy.
The third party software used in this study is dutifully identified in methods section together with the corresponding versions.
The code used to perform background field removal as a part of the phase processing is being prepared to be made publicly available.

Acknowledgments
This study was supported by the German Research Foundation ( DFG, DE2516/1-1 , RE1123/21-1 ) and the Austrian FWF ( I 3001-B27 ). We would like to thank all those working on the variety of open source projects that made this work possible, including the software cited above and the broader scope of the Python scientific stack that did not receive an explicit citation.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.neuroimage.2021.118442 .

A. Mathematical transformation of the orientation dependence expression
The expression for the apparent susceptibility of a white matter fibre, given by W.  , is: where f lip is the volume fraction of myelin lipids, χ lip ∥ represents the magnetic susceptibility of lipid molecules along the direction of the lipid chain, i.e. in the radial direction of the cylinder, while χ lip ⊥ is the susceptibility in the direction perpendicular to the lipid molecule. For convenience, we rewrite Equation A1 with co s 2 θ as: (A2) To simplify the notation, we introduce the magnetic susceptibility components with respect to the fibre orientation (by substituting θ = 0 • and θ = 90 • in Equation A1 ): Note that for brevity, we omit the superscript "fib" in the main text and in Equation A1. We can also relate the mean magnetic susceptibility (MMS) as defined in ( X.  for the susceptibility with respect to the fibre , to the mean magnetic susceptibility of the lipid molecules used in ( W.   (Not surprisingly, as the trace of the susceptibility tensor, MMS remains invariant to the change of the coordinate system and is only scaled by the volume fraction of the lipids.) For completeness, we mention that the magnetic susceptibility anisotropy defined in the present work can be written as follows:

B. Corpus callosum segmentation
1 The probability map of the commissural fibres assigned to CC by TractSeg was thresholded at 90% probability. 2 The mask for the lateral ventricles was extracted from the segmentation by Freesurfer (labels 4 and 43) and dilated with a 6-member structural element. 3 The mask of the central sagittal plane of CC was extracted from the segmentation by Freesurfer (labels 251 to 255) and the mask of the dilated lateral ventricles was removed from it to reduce the effect of possible segmentation errors. The resulting mask was additionally opened with a 6-member structural element. 4 The span of the final mask of CC was defined by morphological dilation of the sagittal plane obtained under item (3) with an ellipsoid with the semiaxes of 11 mm in LR direction, 11 mm in AP direction, and 7 mm in IS direction. 5 The intersection of the dilated sagittal plane and the commissural fibre mask was further morphologically opened with a 6-member element and the largest contiguous component of the resulting mask was taken. 6 Subdivision into the genu, body, and splenium was based on expanding the labels 255, 252 to 254, and 251 correspondingly from Freesurfer along the LR axis

C. Choroid plexus segmentation
1 The mask of the lateral ventricles including the choroid plexus was obtained with Freesurfer. 2 R * 2 values were clipped to the limits of [ 0 , 20 ] s −1 within the mask and set to 0 elsewhere. 3 The R * 2 maps were median-filtered with a spherical kernel of radius 1 mm. 4 Otsu's method ( Otsu 1979 ; as implemented in scikit-image) was applied to identify the optimal threshold to segment out the choroid plexus. 5 Voxels with R * 2 greater than the threshold were considered to be part of the choroid plexus. 6 The resulting mask was further opened with a spherical kernel of radius 1 mm and removed from the mask of the lateral ventricles.