Static and dynamic BOLD fMRI components along white matter fibre tracts and their dependence on the orientation of the local diffusion tensor axis relative to the B 0 -field

Baseline components

The T₂* value of myelinated tissue depends on the orientation of the myelin relative to the B₀ axis.^25–28 Similarly, the T₂* value of a voxel containing tissue and a vessel filled with blood depends on the orientation of the vessel relative to B₀. While the exact effect depends on several factors—such as whether the myelinated structure or blood vessel is centred within the voxel, or perhaps outside of the voxel, the partial-volume/blood-volume fraction, the true shape of the voxel including its point-spread function and aspect ratio, and other factors ^{29 ,30}—the changes in T₂* value with orientation with all other factors held constant takes on a simple form. If we represent the angle between the axes of the vessel and the B₀ field as θ, where θ represents the polar angle in a spherical coordinate system centred around the B₀ axis, the baseline T₂*-weighted signal level of the voxel containing the vessel, denoted S(θ), exhibits an orientation dependence that follows a cos ² θ relationship. ³¹ Most conventional fMRI data are physiological-noise-dominated ³² and for these data the T₂*-weighted signal changes, σ_b, scale positively with the baseline, i.e., σ_b = λ ⋅ S, where λ is a proportionality factor, as is well-known from classic studies of physiological noise^33–35; if this assumption holds, then σ_b ∝ λ ⋅ cos ² θ + β, where β is an offset, and so both the T₂*-weighted baseline signal S and the signal change σ_b will both be higher for voxels that contain myelinated structures and/or vessels parallel to B₀, and both will be lower for voxels that contain myelinated structures and/or vessels that are perpendicular to B₀ as a consequence.

Fluctuation components

Intra- and extravascular field variations induced by changes in blood flow, volume or deoxy-haemoglobin concentration also depend on the orientation of the vessel relative to B₀. Extravascular effects are thought to dominate BOLD for field strengths ≥3T ¹⁹ ; a vessel that is aligned perpendicular to B₀ will exhibit a maximum extravascular field perturbation in its vicinity for a given change in blood oxygenation, whereas a vessel that is aligned parallel to B₀ will induce no extravascular field offset.^16,36 Hence, the dynamic fluctuations of the BOLD signal will also exhibit an orientation effect. We define σ_f as the fluctuation amplitude due to dynamics, which follows a sin ² θ relationship^20,21 and therefore scales in a way that is complementary to the baseline component σ_b defined above. This dynamic fluctuation effect is directly proportional to the orientation angle as σ_f ∝ sin ² θ and thus σ_f vanishes for voxels containing vessels parallel to B₀ with θ = 0°.

MRI data

Data were acquired using protocols approved by the Washington University institutional review board. Written informed consent was obtained from all subjects by the HCP according to the Declaration of Helsinki. A total of 110 subjects from the HCP 3T dataset (S1200, release 03/2017) were chosen at random, avoiding subjects flagged by the HCP for quality control issues.^37–39 These data were acquired on a 3T system based on the Skyra (Siemens Healthineers, Erlangen, Germany) with a 32-channel head coil array; details have been published elsewhere. ³⁹ The “unprocessed”, anonymised image data were downloaded. ⁴⁰

Functional data

The rs-BOLD-fMRI data consisted of four scans (two sets of opposite phase-encode directions, right-to-left and left-to-right) of gradient-echo simultaneous multi-slice (SMS) EPI with blipped CAIPI^41–43 with 2-mm isotropic nominal resolution, matrix = 104 × 90, 72 slices, multi-band factor (MB)=8, no partial Fourier, no in-plane acceleration, TR = 720 ms, TE = 33.1 ms, echo spacing = 0.58 ms, readout bandwidth (BW)=2290 Hz/pixel, time points per run = 1200, total acquisition time ∼15 minutes. One single-band (i.e., MB = 1) volume was acquired in each run as a high-tissue-contrast reference for motion and distortion correction and registration to the T₁-weighted structural data.

Diffusion data

Diffusion-weighted MRI data consisted of six runs of spin-echo SMS EPI. These dMRI data were acquired with three distinct sets of diffusion directions, each acquired with opposite phase-encode directions, right-to-left and left-to-right, with 1.25-mm isotropic nominal resolution, matrix = 168 × 144, 111 slices, MB = 3, partial Fourier = 6/8, TR = 5520 ms, TE = 89.5 ms, echo spacing = 0.78 ms, BW = 1488 Hz/pixel, FA = 78°; refocusing FA = 160°, b-values = [1000, 2000, 3000] s/mm², total acquisition time ∼55 minutes.

Data pre-processing

Structural data pre-processing

The T₁-weighted anatomical data were processed with the minimal pre-processing of the HCP (pipeline version 3.25.0), ⁴⁴ which includes corrections for gradient non-linearities and intensity bias-fields, followed by automatic brain parcellation using FreeSurfer ⁴⁵ version 5.2 which provided the WM segmentation used to define the WM mask in our fMRI and dMRI data.

Functional data pre-processing

The HCP fMRI pre-processing consists of temporal domain operations, such as detrending, and spatial operations, such as motion estimation and correction, geometric distortion correction (for both susceptibility-induced EPI distortion and gradient-nonlinearity distortion), intensity bias correction due to transmit/receive coil nonuniformities, and spatial normalization operations such as registration to atlas space. The pipeline combines all spatial transformations to perform a single one-step resampling to minimize interpolation errors. ⁴⁶ We adapted the pipeline to omit the nonlinear spatial normalization and retain the fMRI data in subject-native space, because this operation rotates the head position and can change the local brain orientation, which complicates proper calculation of orientation at the time of fMRI acquisition. Rigid registration between the T₁-w anatomical data and the distortion-corrected EPI reference scans was calculated with FLIRT from the FSL software package (https://fsl.fmrib.ox.ac.uk). ⁴⁷ We used these transformations to align the DTI tensor fit from its native space into the native functional space corresponding to the head position and pose during each fMRI run to calculate orientation angles during the BOLD-fMRI session (see below).

After pre-processing, we calculated the temporal mean μ and temporal standard deviation σ for each fMRI voxel. No physiological noise removal was performed.

Fibre orientation estimation

We estimated fibre orientations as the angle between B₀ and the primary eigenvectors of the diffusion tensor as θ = acos ⁡ ( e → B 0 · e → DTI ) . The diffusion tensor model was fitted to the b = 0 and b = 1,000 s/mm² data using the FSL “dtifit” function. To account for crossing fibres, we further downloaded the “bedpostX” data calculated by the HCP. This data estimates fibre orientations using the “ball-and-sticks” model (with three sticks) using the “bedpostX” function of FSL^48,49, which used all b-values. The results provided three fibre orientations and their associated volume fractions for every voxel, with one, two or three fibres estimated. Fibres with volume fractions lower than 5% were excluded.

We then resampled the image volumes containing the fibre orientations in the diffusion image space to the rs-fMRI space using nearest-neighbour interpolation. As noted above, to account for rotations of the head between the dMRI and fMRI scans, and ensure that the orientation dependencies of the fMRI data were assessed using the orientation of each tract at the time of the fMRI data acquisition, the fibre orientation(s) at each voxel were rotated using the registration between the dMRI and fMRI data with the tool “fdt_rotate_bvecs” from FSL. We refer to the DTI-based fibre orientation estimate as θ_DTI and bedpostX-based fibre orientation estimates as θ_bX1 (primary fibre orientation), θ_bX2 (secondary fibre orientation) and θ_bX3 (tertiary fibre orientation) rotated into the native coordinate frame of the acquired fMRI data. All angles θ represent orientation as the polar angle in spherical coordinates with the scanner B₀ axis as the zenith direction, which yields a more straightforward interpretation of orientation effects compared to, e.g., representing orientation using Euler angles. ²³

WM and fibre bundle mask generation

For each subject we first created a binary WM mask in native space from the FreeSurfer WM parcellation derived from the MPRAGE anatomical data. We eroded the mask with two iterations of the default 3 × 3 × 3 kernel two remove an approximately two-voxel layer of boundary voxels from the mask to avoid partial volume effects with the cortical GM. To generate individual fibre bundle masks, we used the labelled streamline-density atlas of the IIT Human Brain Atlas. ⁵⁰ This atlas contains 24 tract labels. We registered the atlas reference image to each subject’s fMRI reference image via the subject’s T₁-weighted anatomical reference to map the atlas into subject native space. To create tract-specific masks we then binarized the streamline-density map at a low threshold of 3% of the maximum streamline density value for each tract, as recommended ⁵⁰ ; note that tracts can overlap, i.e., these masks are not mutually exclusive. These masks were used for regional analysis (described below) of the fMRI data based on the DTI-based fibre orientation estimates θ_DTI or bedpostX-based fibre orientation estimates θ_bX1, θ_bX2 and θ_bX3.

Characterization of orientation effects

Our goal was to investigate the dependence of the baseline signal μ and the signal fluctuation σ on fibre orientation to B₀. For each subject and within each bundle mask, we sorted μ and σ values by fibre orientation angle θ in bins of 5° width. Effects are expected to be maximised and minimised at angles 90° and 0°, respectively. Thus we normalised each distribution to the value of the bin closest to 90° (i.e., the bin at θ = 87.5 ± 2.5°). We refer to these distributions as orientation profiles μ(θ) and σ(θ). For each subject we averaged the orientation profiles of all four rs-fMRI runs. Subsequently we averaged orientation profiles over all 110 subjects. We computed these orientation profiles separately for all fibre orientation estimates θ_DTI, θ_bX1, θ_bX2 and θ_bX3 using their respective masks. Histograms of the orientations contained within each fibre mask were calculated, and are presented as Supplemental Material (Supp. Fig. S1). From the initial 24 bundle masks, 16 exhibited sufficient ranges of orientations and were retained for further analysis; the remaining 8 masks, for example the fornix, yielded orientation histograms with too few counts in several bins to adequately cover the full range from parallel to perpendicular orientations relative to B₀. Fibre mask sizes ranged between an average of 801 and 2697 voxels contained within each bundle; a list of the voxel counts for all bundle masks included in our analysis is provided as Supplemental Material (Suppl. Table S1).

Within each retained fibre bundle mask, we calculated the maximum orientation bias between voxels with perpendicular and parallel orientations as Δμ = μ _{θ =90°} − μ _{θ =0°} for the signal baseline; and Δσ = σ _{θ =90°} − σ _{θ =0°} for the fluctuation amplitudes, respectively. We refer to these quantities as signed orientation signatures. Signal baseline components σ_b are expected to exhibit a signed orientation signature greater than zero since σ_b ∝ cos ² θ, thus we refer to masks with a positive trend Δσ > 0 as “baseline-dominated”, whereas signal fluctuation components σ_f are expected to exhibit a signed orientation signature less than zero since σ_f ∝ sin ² θ, thus we refer to masks with a negative trend Δσ < 0 as “fluctuation-dominated”. We performed this analysis using all voxels within each WM mask, and then, to test whether a stronger effect may be seen in voxels with more consistent fibre orientations, we repeated the analysis on only those voxels whose fractional anisotropy (FA) exceeded a moderate threshold (FA ≥ 0.2) and high threshold (FA ≥ 0.5) threshold. We focus however on the results from the full WM masks without FA thresholding, since these are expected to be more conservative, with potentially stronger orientation effects expected when considering higher FA thresholds.

Orientation effects based on fibre direction estimates from DTI fits

Orientation effects from pooling across the entire WM

The baseline signal μ, averaged over the entire WM, exhibits a clear dependence on the local WM orientation estimated from the dMRI data. Voxels for which the dominant diffusion direction was parallel to B₀ exhibited an ∼11% higher baseline T₂*-weighted signal than voxels for which the dominant diffusion direction was perpendicular to B₀, as shown in Figure 3(a). This baseline signal change corresponds to a change in transverse relaxation rate, ΔR₂*, of 3.1 Hz between voxels in which the principal direction is parallel to or perpendicular to B₀ assuming a mono-exponential signal decay, and for our TE of 33.1 ms. The signal fluctuation σ, however, averaged over the entire WM shows a mixed bias relationship, as shown in Figure 3(b), where error bars are of comparable magnitude to the effect itself.

OPEN IN VIEWER

Figure 3.

The orientation dependence of the fMRI time series within the whole WM mask and averaged over 110 subjects. Error bars are the standard deviation across subjects. Here the fibre directions are estimated from DTI fits to the dMRI data (and thus are denoted as θ_DTI). (a) The baseline signal as a function of the orientation of the fibre direction estimate relative to the B₀ axis shows an orientation bias of about 11% between voxels that contain perpendicular and parallel fibres and (b) The fluctuation amplitude of the time series shows a bi-phasic behaviour and exhibits error bars of similar magnitude as the effect, suggesting no effect or a spatially heterogenous distribution of orientation effects.

Orientation effects observed within individual fibre bundle masks

Heterogenous orientation profiles emerge when the dynamic BOLD signal changes are plotted for individual fibre bundles. Profiles are displayed in Figures 4 and 5. All fibre bundle masks tested exhibited clear orientation effects in their baseline T₂* signal μ, as expected, which matches our observations from the whole WM mask. The consistency of this observed effect suggests that all fibre bundle masks tested contained a sufficient number and range of orientations to enable detection of orientation effects in our T₂*-weighted data.

OPEN IN VIEWER

Figure 4.

Eight (out of 16) fibre bundles exhibit negative orientation bias profiles as expected from fluctuation changes in BOLD. The profiles are averaged over 110 subjects and error bars represent the standard deviation across subjects. Orientation angles refer to the principal diffusion direction estimate from the DTI tensor fit θ_DTI. Blue plots are the orientation profile of the times series temporal mean μ (baseline) and black plots are the temporal standard deviation σ (fluctuation). Fibre bundles are sorted with respect to the amount of negative fluctuation bias Δσ from highest to lowest bias (a–h); exact values are in Table 1.

OPEN IN VIEWER

Figure 5.

Three (out of 16) fibre bundles (a–c) exhibit positive orientation bias profiles that suggest that their time series fluctuations are dominated by baseline fluctuations. Five fibre bundles (d–h) exhibit no orientation bias beyond the magnitude of the error bars. The profiles are averaged over 110 subjects and error bars represent the standard deviation across subjects. Orientation angles are the principal diffusion direction estimate from a DTI tensor fit θ_DTI. Blue plots are the time series temporal mean μ (baseline) and black plots are the temporal standard deviation σ (fluctuation). Fibre bundles are sorted with respect to the amount of positive fluctuation bias Δσ from highest to no bias (a–h); exact values are in Table 1.

Fibre bundles are categorised by their orientation signature in Table 1. Eight of the 16 bundles are classified as fluctuation-dominated with a negative orientation signature Δσ in their fluctuation amplitudes; their corresponding profiles are plotted in Figure 4. Three bundles are classified as baseline-dominated with a positive orientation signature Δσ; their corresponding profile are plotted in Figure 5(a) to (c). Five bundles are classified as bias-free with little-to-no orientation bias beyond the magnitude of the error bars; their corresponding profile are plotted in Figure 5(d) to (h).

OPEN IN VIEWER

Table 1.

Fibre bundles categorised by their orientation bias signature Δσ.

Fibre bundle	Baseline orientation bias Δ μ ( μ θ = 0 ° − μ θ = 90 ° ) [%]	Fluctuation orientation bias Δ σ ( σ θ = 0 ° − σ θ = 90 ° ) [%]
Inferior fronto-occipital fasciculus	18.2	–19.3
Corpus calossum forceps major	12.6	–18.2
Optic radiation	17.5	–17.4
Corpus calossum forceps minor	10.6	–15.6
Inferior longitudinal fasciculus	18.2	–14.8
Corpus calossum middle	6.1	–8.9
Cingulum	6.6	–6.8
Vertical occipital fasciculus	11.0	–6.1
Corticospinal tract	8.8	21.1
Parieto-pontine tract	8.0	20.0
Frontal aslant tract	9.8	9.4
Middle longitudinal fasciculus	8.1	–5.1
Occipito-pontine tract	9.1	–5.0
Fronto-pontine tract	9.6	4.2
Arcuate fasciculus	10.0	–3.8
Superior longitudinal fasciculus	9.2	0.3

Note: Orientation angles in this analysis are the principal diffusion direction estimate from a DTI tensor fit θ_DTI. The first eight bundles exhibit a negative orientation bias Δσ as is expected from fluctuation BOLD fluctuations from voxels containing vessels running parallel to the fibre directions, their profiles are in Figure 4. Three fibres exhibit a positive orientation bias Δσ as is expected in the case of pure baseline dominated signal fluctuations, their profiles are in Figure 5(a) to (c). Five fibres show no clear orientation bias, their profiles are in Figure 5(d) to (h).

The bundles with the strongest negative fluctuation bias Δσ corresponding to the fluctuation-dominated signature, also exhibit the strongest baseline biases Δμ. The most fluctuation-dominated bundles are the inferior fasciculus, the corpus callosum and the optic radiation, with Δσ from −15–19%, and Δμ from +11–18%. The heterogenous profiles between the fluctuation-dominated bundles (Figure 4), baseline-dominated and bias-free bundles (Figure 5) may explain the bi-phasic σ(θ) profile seen in Figure 3(b) obtained when averaging over the entire WM mask. The results did not change in their interpretation when we repeated the analysis with more restrictive WM masks based on selecting only voxels with FA values above either 0.2 or 0.5 (see Supplementary Material, Suppl. Table S2).

Orientation effects based on fibre direction estimates from probabilistic tractography

It is well-known that fibres can cross and kiss within a given voxel, in which case the principal direction of the DTI fit will be an inaccurate estimate of the true fibre directions and thus the inferred vessel directions. We therefore performed the same orientation profile analysis for the individual fibre masks using three direction estimates from probabilistic tractography (Figure 6) which can resolve multiple fibre crossings in a voxel.

OPEN IN VIEWER

Figure 6.

Exemplar orientation profiles from the probabilistic direction estimates from bedpostX. The profiles for the major direction is plotted in dark blue (μ(θ_bX1)) and black (σ(θ_bX1)); the results for the intermediate direction is in light blue (μ(θ_bX2)) and light grey (σ(θ_bX2)); and the results for the minor direction is in lightest blue (μ(θ_bX3)) and lightest grey (σ(θ_bX3)). (a, b) Dynamics-dominated bundles of the corpus callosum and the optic radiation; (c, d) baseline-dominated examples of the parietopontine tract and the frontal aslant tract. Trends are similar to the DTI-based fibre direction angle estimates shown in Figures 4 and 5 when compared with the major direction but diminish for the intermediate and the minor fibre direction estimates.

Orientation effects in the gradient-echo BOLD baseline signal due to baseline T2* effects

We found that the baseline gradient-echo T₂*-weighted signal, calculated as the temporal mean of the T₂*-weighted time-series data, varies on average by 11% between WM voxels that contain parallel and perpendicular fibre direction estimates. Our data consist of T₂*-weighted images, and not quantitative T₂* maps, and while other non-T₂* components may exhibit subtle orientation effects such as any T₁ weighting^51,52, if we assume that the orientation effects are dominated by T₂* our observation corresponds to a change in relaxation rate, ΔR₂*, of ∼3 Hz. Similar observations have been reported based on fibre orientation estimates derived from DTI: a 3T study by Denk et al. reported parallel-to-perpendicular orientation ΔR₂* changes of approximately 2–3 Hz in the WM using multi-echo susceptibility imaging ³¹ and Kaczmarz et al. reported a ΔR₂* of 2.5 Hz in a 3T study. ²³ Increasing field strength will enhance the orientation dependence, and indeed observed ΔR₂* changes of ∼6 Hz were found in the corpus callosum in two ex-vivo studies at 7T.^25,27 These studies also indicate that effects of blood was only a minor contributing source of orientation effects, since these brain-tissue specimens presumably contained far less blood than what is present in vivo. The latter ex-vivo study also found that ΔR₂* did not change notably after tissue iron had been extracted from the specimen, which suggests that the observed dependence may be driven predominantly by myelin. Also a recent study investigating orientation effects in spin-echo imaging demonstrated that observable orientation effects can be generated by myelin alone, ⁵³ suggesting that blood vessel effects are not needed to generate an orientation dependency. Here, we show how the baseline orientation dependence is heterogenous between distinct fibre tracts, with Δμ ranging between 6% to 18% across bundle masks. This can in part be attributed to the heterogeneous myelination levels across the brain, causing some tracts to be more myelinated than others such that they exhibit varying levels of orientation effects. Still, another source that should contribute to the heterogeneity in the baseline signal orientation dependence in our data are blood vessels—in particular veins that contain T₂*-shortening deoxygenated blood—that run parallel to many of the bundles.

Support of a role of the vasculature on the observed orientation effects is provided by a recent study by Hernández-Torres et al. on the orientation dependence of blood volume and flow estimates in the WM derived from DSC imaging. By administering a contrast agent, they found an orientation-dependent change in R₂* that was consistent with approximately half of the blood volume in the WM being within blood vessels that were anisotropic and aligned parallel to the local WM fibre direction; they also found estimated blood volume to be almost 100% higher in perpendicular WM than in parallel WM at the highest contrast agent dose. ²² A similar conclusion was reached by Doucette et al., who, using spin-echo based DSC data, also estimated from the observed orientation effect that roughly half of the WM vascular volume is comprised of vessels running in parallel with WM fibre tracts, ²⁴ which further suggests that blood vessels with smaller diameters seen in spin-echo data also exhibit anisotropy and are aligned with the WM tracts; in their case, numerical simulations indicated that the spin-echo signal was sensitive to small vessels with radii of about 100 µm on average. An orientation effect was also found in a later study by Polimeni et al. who used high-resolution DSC imaging and DTI data at 7T to assess the B₀-orientation effect in the estimation of cerebral blood volume. ⁵⁴ This study confirmed that even when using small voxels the blood volume estimates show an orientation bias: on average, voxels that contained parallel fibres yielded a ∼10% lower blood volume estimate than voxels that contained perpendicular fibres. In agreement with our findings, this recent study also found heterogenous orientation profiles in DSC-based blood volume estimates across fibre bundles using a similar analysis approach, supporting the notion that vascularisation across the WM bundles might be heterogenous. Because partial volume effects between vessels and surrounding tissue within a voxel play an important role on the orientation effect, ³⁰ and interactions between voxel geometry and extravascular fields can impact observed signal levels, ²⁹ especially for voxels containing vessels parallel to B₀, further work will be needed to fully account for the observed differences in the orientation bias across imaging resolutions and field strengths.

Orientation effects in gradient-echo BOLD fluctuation amplitudes and implications for WM-fMRI

The novel aspect of our work was to exploit the orientation dependence of T₂* to test for potential dynamic sources in gradient-echo BOLD amplitude fluctuations that could be more directly associated with blood vessels running alongside WM fibre tracts. Static baseline signal components, from sources such as myelin and venous blood, should introduce an orientation bias Δσ that tracks the baseline orientation bias Δμ that would manifest as a maximal T₂* change for tracts parallel to the B₀ axis. However, if there are dynamic changes in blood oxygenation levels, i.e., BOLD effects, within the vessels that dominate the orientation effect, the orientation bias should follow an opposite negative trend that would manifest as a maximal BOLD fluctuation for tracts perpendicular to B₀. When pooling across the entire WM, we found an orientation effect size that was small relative to the error bars, however the relationship between BOLD fluctuation amplitude varied slowly with orientation, suggesting that there may be a systematic relationship. A more apparent orientation effect was observed by considering individual fibre bundles. Here we found that about half of the fibre bundles show this fluctuation-dominated orientation signature consistent with an effect driven by blood vessels. The strongest effects were seen in the inferior fronto-occipital and longitudinal fasciculi, the corpus callosum and the optic radiation. Given the recent reports of WM-fMRI activity it might be of interest to study these fibre bundles in particular, as BOLD responses from these bundles may be more readily detectable. In fact many WM-fMRI studies have reported BOLD activity in the corpus callosum.^1–3 Although our results are based on observations of gradient-echo BOLD data, it is expected that spin-echo BOLD data would show similar orientation biases based on previous reports, mentioned above, that suggest that blood vessels far smaller than the voxel size also appear to run alongside WM tracts and generate measurable orientation effects. ²⁴

Our study is limited by sensitivity variations in the fMRI data across the cerebrum which could impact the detectability of fMRI fluctuations and could be a reason why we do not detect BOLD fluctuations in some fibre bundles. An additional limiting factor is the high slice-acceleration (MB = 8) used in the fMRI protocol. Such high parallel-imaging acceleration contributes to spatially-varying noise levels, with reductions of SNR in the centre of the head including some WM regions, thus we expect effects to be more detectable in studies using less acceleration. Furthermore, partial volume effects could also hinder our ability to detect BOLD fluctuations in all tracts of the WM. However, our results are in agreement with a recent contrast-enhanced vascular imaging study by Bernier et al. who ranked fibre bundles by vascularisation using the same fibre bundle atlas as in this work. ⁵⁵ The bundles that appeared most vascularised in their data are also the bundles that we find to have the highest dynamic BOLD fluctuations. This does not imply that the bundles where neither vessels (in their data) nor BOLD fluctuations (in our data) were detected are not vascularised per se, but it may mean that these particular bundles have reduced blood volume, making detection both in this previous study and in our study more challenging.

Crossing fibres are a confounding factor to estimate fibre orientations from dMRI data and thus to infer blood vessel orientation. For example, the corpus callosum is more uni-directional than other fibre bundles, which potentially renders it more favourable for accurate estimation of fibre direction and thus more accurate blood vessel orientation estimation, which also could improve the detectability of BOLD orientation effects there. However, our results remained mostly consistent when we used the multi-fibre estimates obtained with a crossing-fibre model (bedpostX) instead of the single-fibre orientation from the diffusion tensor model. Our multi-fibre analysis led to the bias signature of only two fibre bundles, namely the parietopontine and the occipitopontine tracts, changing between the primary and tertiary orientation estimates. A potential explanation is that the fibres corresponding to the primary orientation possess no parallel running vessels, whereas the fibres along the tertiary orientation do possess parallel vessels that exhibit a dynamic BOLD component.

WM metabolism and sources of dynamic BOLD

Our results suggest that parallel-running vessels in certain bundles carry dynamic BOLD fluctuations in the resting state. The question still remains regarding which processes drive the changes in blood oxygenation levels in the WM vasculature. One possibility is that these fluctuations reflect global effects such as physiological noises in WM driven by the respiratory/cardiac cycles. Previous fMRI studies presented evidence for physiological noise within WM,^8,56 however these fluctuations are expected to be global as they are driven systemically, yet in our data fluctuations are only observed in a subset of WM tracts. Another possibility is that the observed BOLD fluctuations reflect local effects such as functional hyperaemia driven by local fluctuations in metabolic demand. Unfortunately, knowledge about WM metabolism and neurovascular coupling is sparse. GM metabolism is well understood, and it is well accepted that oxidative phosphorylation is the dominant energy source. GM energy is suggested to be mainly used for neural signalling compared to non-signalling processes. ⁵⁷ In contrast, WM energy is mainly consumed for non-signalling “housekeeping” purposes, such as maintaining homeostasis, and overall the WM energy consumption is estimated to be three times smaller than that of the GM, ⁵⁸ which is roughly consistent with the observation that the WM exhibits a vascular density that is four times lower than that of GM. However, while less than half of energy the consumption in WM is related to signalling, this is still relatively “high” given that synapse density is substantially lower; to our knowledge this has only been measured in rat brain, where synapse density is 80× lower in WM compared to GM ⁵⁸ , and collectively WM synapses require less than 0.5% of the energy of GM synapses. ⁵⁸ Rather than synaptic activity, other signalling processes such as the propagation of action potentials along axons and the maintenance of their resting membrane potentials appear to drive metabolic demand in WM. Still, there is no consensus on how much metabolism in the WM is based purely on aerobic glycolysis, i.e., glycolysis without oxidative metabolism when there is no lack of available oxygen. Animal studies are somewhat contradictory, with one rodent study reporting higher glycolytic than oxidative metabolism in WM, ⁵⁹ another study reporting a correlation between WM oxidative glucose metabolism and GM electrostimulation ⁶⁰ and a more recent study suggesting that axonal oxidative glucose metabolism as well as lactate production (a by-product of glycolysis) rise when activity is evoked in the GM regions sending or receiving the associated WM projections. ⁶¹ In summary, it might be that WM utilises a mix of metabolic mechanisms, of which some require oxygen supply via the vasculature and thus requiring functional hyperaemia and the consequent BOLD effect. Whether BOLD fluctuations seen in the WM are truly linked to neuronal signalling in GM and subsequent activation of WM axons cannot be answered here, and the BOLD effects we have observed could be entirely unrelated to neural activity.

Physiological noise removal was not performed here in part because our hypothesis was that BOLD fluctuations in the WM should exhibit orientation effects, and these fluctuations can be driven by neuronal activity or by systemic physiology. We previously demonstrated that some physiological noise components that reflect true oxygenation changes exhibited orientation dependence in the cortical GM, whereas other physiological noise components that reflect other non-BOLD changes did not exhibit orientation dependence in the cortical GM, as expected. ²¹ Neuronally-driven fluctuations in the BOLD signal exhibit an orientation dependence, and so while physiological noise correction, if performed properly, will have some impact on the size of the orientation dependence effect. We note that our derivation for why baseline-dominated orientation dependence given by σ = λ ⋅ S ∝ cos ² θ assumes that the fMRI data in the physiological noise dominated regime, and so fMRI data with lower levels of physiological noise are expected to more often exhibit a fluctuation-dominated orientation dependence given by σ ∝ sin ² θ. It is also unclear how “draining vein” effects, which cause a displacement of BOLD signal changes away from the site of neural activity, manifest in the WM. Although the orientation dependency is expected to be dominated by larger (i.e., draining) veins—so by definition there will be some influence of draining veins—proper draining-vein displacement effects could be observed far downstream from the capillary bed, e.g., from penetrators draining subcortical GM structures, or an increased orientation effect at the downstream end of a tract caused by an accumulation of deoxygenated blood from all upstream capillary beds. Note that most of the venous blood of the WM is drained outwards, and so deoxygenated blood from the WM may be observed downstream in the cortical GM. Still, observed fMRI fluctuations along tracts might offer a means to study vascular function and integrity in the WM. This might be relevant for research into neurodegenerative diseases, such as vascular dementia where WM tissue pathology and vascular impairment occur as concomitant phenomena.^62–64