Spatiotemporal topological correspondence between blood oxygenation and glucose metabolism revealed by simultaneous fPET-fMRI in brain’s white matter

Abstract

White matter (WM) makes up half of the human brain. Compelling functional MRI evidence indicates that white matter exhibits neural activation and synchronization via a hemodynamic window. However, the neurometabolic underpinnings of white matter temporal synchronization and spatial topology remain unknown. By leveraging concurrent [¹⁸F]FDG-fPET and blood-oxygenation-level-dependent-fMRI, we demonstrated the temporal and spatial correspondences between blood oxygenation and glucose metabolism in the human brain white matter. In the temporal scale, we found that blood-oxygenation-level-dependent signals shared mutual information with FDG signals in the default-mode, visual, and sensorimotor-auditory networks. For spatial distribution, the blood-oxygenation-level-dependent functional networks in white matter were accompanied by substantial correspondence of FDG functional connectivity at different topological scales, including degree centrality and global gradients. Furthermore, the content of blood-oxygenation-level-dependent fluctuations in the white matter default-mode network was aligned and liberal with the FDG graph, suggesting the freedom of default-mode network neuro-dynamics, but the constraint by metabolic dynamics. Moreover, the dissociation of the functional gradient between blood-oxygenation-level-dependent and FDG connectivity specific to the white matter default-mode network revealed functional heterogeneities. Together, the results showed that brain energy metabolism was closely coupled with blood oxygenation in white matter. Comprehensive and complementary information from fMRI and fPET might therefore help decode brain white matter functions.

Introduction

Since its inception, blood-oxygenation-level-dependent (BOLD)-fMRI has generated interest as a tool for detecting functional activity in cerebral gray matter (GM), and as a means of studying the neuro-dynamics in white matter (WM) (Gawryluk et al. 2014). WM is a long-neglected half of the brain (Fields 2013), with its BOLD signals regarded as “noise” and deliberately regressed out (Caballero-Gaudes and Reynolds 2017; Grajauskas et al. 2019). However, compelling evidence has shown that changes in WM BOLD signals may reflect neural activities (Ji et al. 2017; Gore et al. 2019). Consequently, there has been an increased interest in the neurobiology of functional fluctuations within the WM. This has ranged from neuroanatomy (Ji et al. 2017) to studies of neuropathology (Ji et al. 2019; Fan et al. 2020; Li et al. 2020b), neurobehavior (Li et al. 2020a), neuro-topology (Peer et al. 2017; Li et al. 2019), neurophysiology (Greene et al. 2020; Li et al. 2021b), and neurogenetics (Li et al. 2021d). Increased knowledge of the functional role of WM has been investigated by BOLD-fMRI that reflects the complex interactions of blood oxygenation and changes in blood flow, blood volume, or intravascular magnetic susceptibility (Tsvetanov et al. 2021). Knowledge concerning the associations between hemodynamics and underlying metabolic dynamics is less understood.

Characterization of the glucose metabolic basis of BOLD has been central to understanding WM's energetic dynamics across multiple levels, including temporal synchronization and spatial correspondence (Harris and Attwell 2012; Gawryluk et al. 2014). To this end, PET shows activity-dependent glucose metabolic changes in the WM, where glucose is the primary energy source (Yu et al. 2018). PET measures an intermediate stage in the chain linking neural activity via metabolism to the BOLD mechanism and provides complementary information because of its high sensitivity and ability to track specific biochemical processes (Frackowiak et al. 2004). Using the PET tracer [¹⁸F]-fluorodeoxyglucose ([¹⁸F]FDG), it was initially assumed that FDG uptake measured neuronal glucose consumption (Li et al. 2021a). However, recent studies have reported that activities of glial cells (high glia-to-neuron ratio in the WM), including astrocytes (Zimmer et al. 2017; Carter et al. 2019; Bonvento and Bolanos 2021; Howarth et al. 2021), microgliacytes (Xiang et al. 2021; Zimmer et al. 2022), and oligodendrocytes (Harris and Attwell 2012; Freeman and Rowitch 2013; Howarth et al. 2021), also consumed substantial glucose. Thus, the primary use of energy in the WM is not for the efficient propagation of action potentials but rather for maintaining the proper functions of glial cells (Gore et al. 2019). Using the [¹⁸F]FDG-PET technique, clinical studies have found that glucose metabolism in the WM was related to some neuropathologies (Buchsbaum et al. 2007; Pagani et al. 2014; Chassoux et al. 2017; Mitelman et al. 2018). Because hemodynamics and glucose metabolism in the WM may relate to neuropathologies, understanding their relationships could provide new phenomenological models of brain disorders.

The feasibility of simultaneously using [¹⁸F]FDG-fPET/fMRI of the human brain has recently been reported to comprehensively decode brain functions and networks (Wehrl et al. 2013; Riedl et al. 2016; Savio et al. 2017; Shokri-Kojori et al. 2019). However, most of these studies acquired static snapshots of brain activations. But fPET is a nascent technique for tracking continuous glucose uptake, which has recently shown promise for assessing the dynamics of neural metabolism (Ionescu et al. 2021; Jamadar et al. 2021). This methodology has provided direct evidence that BOLD activity in the WM was directly related to variations in glucose metabolic demands (Guo et al. 2022). However, this study focused on local activity along with connectivity within the WM, so examining brain WM as an integrative network of functionally interacting regions might provide new insights into larger-scale neural communication (Palombit et al. 2022).

Leveraging available simultaneous [¹⁸F]FDG-fPET/fMRI data at rest (Jamadar et al. 2020), we investigated whether WM functional connectivity (WMFC) derived from blood oxygenation showed underlying glucose metabolic dynamics. First, we determined whether the spatial organization of BOLD functional connectivity corresponded to that of FDG metabolic temporal correlations from multiple topological scales. The second goal was to determine if we could identify how FDG metabolic connectivity might constrain BOLD fluctuations in the WM. We used a graph signal processing (GSP; Huang et al. 2016) to investigate where and to what extent BOLD fluctuations across the WM were organized in a manner aligned with FDG metabolic connectivity networks. Finally, based on our prior study showing the nonuniform distribution of BOLD functional networks in WM (Li et al. 2019, 2021d), we determined if the core functional gradient of BOLD connectivity was similar to FDG temporal connectivity. Gradients constitute a principal organizing axis of the human brain (Margulies et al. 2016; Huntenburg et al. 2018). Studying these correlations with intrinsic dimensions could increase our understanding of how WMFC derived from BOLD relates to glucose metabolic connectivity.

Materials and methods

Study design

The experimental paradigm is shown in Fig. 1A. We analyzed available simultaneous [¹⁸F]FDG-fPET/BOLD-fMRI data from healthy individuals at rest given a continuous infusion of FDG (Jamadar et al. 2020). We first concatenated 6 BOLD sessions into 1 to match the number of segments of FDG. Analogous to the previous approach in constructing BOLD and FDG temporal correlations (Jamadar et al. 2021; Voigt et al. 2022), we focused on voxels within the WM for each subject. We then obtained group-averaged WMFC from BOLD and FDG. It would also be helpful to compare BOLD and FDG WMFC in parallel better to understand the underlying glucose metabolic dynamics of neuro-dynamics. Thus, the main results were focused on temporal synchronization, and spatially topological correspondence based on network perspectives ranging from local nodes to global gradients.

Data acquisition and preprocessing

We used unique simultaneous [¹⁸F]FDG-fPET and BOLD-fMRI data openly available on OpenNeuro (accession number: ds002898). The detailed data are provided by Jamadar et al. (2020; see Online-Only Data Supplement, Supplemental Methods). This design was reviewed by the Monash University Human Research Ethics Committee. The study was carried out according to the guidelines of the Helsinki Declaration of 1975. All subjects provided informed consent to participate in the study.

Each participant underwent a 95-min simultaneous MR-PET scan using a Siemens scanner. The detailed data information and preprocessing are provided in Online-Only Data Supplement, Supplemental Methods.

Creating group-level WM mask

After obtaining each subject's thresholded WM mask from preprocessed data, we identified voxels as WM across all subjects to create a group-level WM mask. To further exclude the impact of deep brain structures, we used a probability (25% threshold) Harvard-Oxford Atlas to remove subcortical nuclei (i.e. bilateral thalamus, putamen, caudate, pallidum, and nucleus accumbens) from the group-level WM mask (Li et al. 2019, 2020a, 2020b, 2021d). Consequently, 12,224 voxels remained in this mask. Subsequent analyses were performed in the abovementioned group-level mask unless otherwise noted.

WM functional systems in BOLD

To determine the WM functional organization of BOLD, we first performed a group spatial independent components analysis (sICA) using GIFT software (version 3.0c, http://icatb.sourceforge.net/; Calhoun et al. 2001). All subjects’ BOLD data were concatenated, and the temporal dimension of the aggregate data set was reduced using principal component analyses within the group-level WM mask (a total of 12,224 voxels). Next, we estimated the time series and spatial maps using an infomax ICA algorithm. Basile et al. (2022) used sICA to identify spatially independent components (IC) in WM based on track-weighted dynamic functional connectivity. They performed sICA decomposition at 3 different dimensional levels by selecting a different number of ICs (n): n = 10 (IC₁₀) for large-scale networks, n = 20 (IC₂₀) as empirical resting-state networks, and n = 100 (IC₁₀₀) for fine-grained parcellation. Although all these ICs showed high intra-subject, inter-subject, and inter-cohort reproducibilities, the IC₁₀ showed the highest test–retest reproducibility. Thus, we used 10 ICs to identify WM functional networks in the present study. The group time series was calculated as the average of the z-scored subject time series. The dual regression (spatial-temporal regression) was used to generate the subject-specific z-score maps (Filippini et al. 2009). The WM functional networks (systems) were identified with a z-score of 1.

Mutual information between BOLD and FDG temporal signals

We next sought to investigate the temporal synchronization between BOLD and FDG signals. The BOLD and FDG temporal signals have different sample points and temporal scales. We thus downsampled the BOLD signals to FDG data. In addition, the FDG reconstruction software could not provide non-integer second fPET bins, resulting in a phase shift in the temporal synchrony of the fMRI and fPET time series (Jamadar et al. 2021). We thus used mutual information (MI) derived from information theory, which quantified their correspondence in WM functional networks for each individual. MI derived from information theory quantified how knowledge of a variable or set of interacting variables reduced the uncertainty of an observed system (Cover and Thomas 1991). This method has been widely used in neural-level (Varley et al. 2023) and circuit-level information processing (Sundaram et al. 2020). Given 2 variables X (i.e. BOLD) and Y (i.e. FDG) from 1 brain voxel, we calculated the MI between the 2 as the degree to which knowing the state of one variable reduced our uncertainty about the state of the other. Formally,

where the entropy, H(BOLD), quantified the randomness of the BOLD signal, and the conditional entropy, H(BOLD|FDG), quantified the randomness of BOLD conditioned on observations of FDG signals.

Construction of BOLD and FDG WMFC

Analogous to a previous approach to constructing BOLD (Liao et al. 2018) and FDG temporal correlations (Jamadar et al. 2021; Voigt et al. 2022) in GM, we only focused on voxels within the group-level WM mask across all subjects. Both BOLD and FDG time series were extracted for each voxel within the group-level WM mask. To construct WMFC, we computed connectivity strength between the time series from any pairs of voxels using Pearson's correlation coefficients (r). Subsequently, a weighted matrix (12,224 × 12,224 voxels) was obtained for each subject and each modality. Finally, we averaged WMFC across all subjects as group BOLD WMFC, as well as FDG WMFC for further assessments (Jamadar et al. 2021).

Degree centrality of BOLD and FDG WMFC

High-degree hubs play a crucial role in integrating information; we therefore first tested the hypothesis that BOLD and FDG WMFC showed a high degree of spatial similarity in degree centrality (DC). The DC is an example of a strictly local measure; it can characterize a single node or voxel (Betzel and Bassett 2017). To minimize the influence of thresholds on DC, we used a series of thresholds (r ranging from 0 to 0.6, and interval = 0.1) of the weighted WMFC matrix. In a given threshold at r_thre, elements within BOLD and FDG WMFC were set to 0 if this element was lower than r_thre (Li et al. 2019, 2020a, 2020b). Then, the still-existing connections for each voxel were summed as the voxel’s DC. Finally, the DC values were obtained for each voxel per modality. Considering the outliers in DC patterns, Spearman's rank correlation coefficient (Rho) was used to investigate the DC spatial similarity between BOLD and FDG.

BOLD fluctuation decompositions into FDG WMFC

Metabolic connectivity has emerged as a complementary approach to decode brain connectivity (Di et al. 2012; Savio et al. 2017; Amend et al. 2019; Jamadar et al. 2019, 2021, 2022). Since FDG-PET is excellent for elucidating the glucose metabolic basis of the BOLD-fMRI signals (Wehrl et al. 2013), we used the GSP method (Huang et al. 2016) to determine the extent of BOLD activity coupling with FDG WMFC (Fig. S1). Brain connectivity is commonly conceptualized as a graph or a network. We first decomposed the functional harmonics of FDG WMFC. The FDG WMFC was defined as the pair |$g=\left(v,\mathrm{W}\right)$|⁠, where |$v=\left\{1,2,\dots, n\right\}$| is a set of |$n$| voxels and |$\mathrm{W}\in{R}^{n\times n}$| represents weights of edges in the network with |${w}_{ij}={w}_{ji}$| being the weight of the edge |$\left(i,j\right)$|⁠, in which |$i,j\in v$|⁠. The Laplacian matrix was defined as the difference |$\mathrm{L}:= \mathrm{D}-\mathrm{W}\in{R}^{n\times n}$|⁠, where |$\mathrm{D}\in{R}^{n\times n}$| was a diagonal matrix, |${\mathrm{D}}_{ii}={\sum}_{j=1}^n{w}_{ij}$|⁠. Given that the graph Laplacian |$\mathrm{L}$| is real symmetric, it can be decomposed into its eigenvalue components |$\mathrm{L}=\mathrm{U}\Lambda{\mathrm{U}}^{\mathrm{T}}$|⁠; the diagonal eigenvalue matrix is defined as |$\Lambda := \operatorname{diag}\left({\mathrm{\lambda}}_0,{\mathrm{\lambda}}_1,\dots, {\mathrm{\lambda}}_{n-1}\right),\mathrm{where}\ {\mathrm{\lambda}}_k$| is the set of eigenvalues and is ordered so that |$0={\mathrm{\lambda}}_0\le{\mathrm{\lambda}}_1\le \dots \le{\mathrm{\lambda}}_{n-1}$|⁠, and |$\mathrm{U}:= \left[{u}_0,{u}_1,\dots, {u}_{n-1}\right]$| is the eigenvector matrix. Spatial patterns of these eigenvectors were regarded as harmonics of FDG WMFC (Glomb et al. 2021). The actual harmonic values are not meaningful per se, but the difference between the values assigned to 2 voxels reflects how different they are in terms of their “FC profile” (Atasoy et al. 2016; Huang et al. 2016; Glomb et al. 2021).

After obtaining functional harmonics of FDG WMFC, we used graph signal filtering to decompose the WM BOLD fluctuations into aligned components (i.e. low-frequency eigenmode) and liberal components (i.e. by higher frequency eigenmode) following previous studies (Atasoy et al. 2016; Huang et al. 2016; Medaglia et al. 2018; Glomb et al. 2021). Aligned components represented signals that varied smoothly across the network (graph), whereas liberal components denoted signals that varied highly across the graph at single moments in time (Medaglia et al. 2018). Because the selection of cutoff frequency methods varied, the spectrum was divided into 3 portions (cutoff frequency was |${C}_1$| and |${C}_2$|⁠) with equal energy based on average energy spectral density (Preti and Van De Ville 2019). The |$n\times n$| matrix |${u}^{\left(\mathrm{low}\right)}$| contained the first |${C}_1$| eigenmode complemented with |$n-{C}_1$| zeros columns, and the matrix |${u}^{\left(\mathrm{high}\right)}$| contained the first |${C}_2$| zeros columns and complemented with |$n-{C}_2$| last eigenmode. Thus, WM BOLD fluctuations were filtered by the following equations (Preti and Van De Ville 2019):

For each subject, the norms of |${\mathrm{s}}^{\mathrm{Aligned}}$| and |${\mathrm{s}}^{\mathrm{Liberal}}$| across times were defined as the aligned and liberal components.

BOLD and FDG WMFC gradient analyses

Connectivity-based gradients capture a global topography (Huntenburg et al. 2018), providing a complementary view to depict brain function at a systematic scale. Rather than marking discrete boundaries at locations of abrupt change in WMFC, we continuously investigated spatial variations across the brain WM along overlapping organizing axes. To characterize the global topography of BOLD and FDG WMFC, we identified spatially gradient axes of inter-voxel functional variations using the BrainSpace toolbox (Vos de Wael et al. 2020; https://brainspace.readthedocs.io/en/latest/pages/install.html). This algorithm starts with computing an affinity matrix between each paired WM voxel. The elements in an affinity matrix represent the degree of similar WMFC of 2 voxels. According to previous studies (Margulies et al. 2016; Dong et al. 2021), the WMFC should be row-wise thresholded before obtaining an affinity matrix, which can affect the similarity calculation. Consistent with our previous studies (Meng et al. 2021; Yang et al. 2021), only the top 10% connections of each voxel remained. Then a cosine distance was computed between any 2 rows of the sparse matrix to obtain a symmetrical similarity matrix. Diffusion map embedding was implemented on the similarity matrix to derive gradients of WMFC. This nonlinear reduction method allowed the connectivity matrix to project into an embedding space. Two parameters control the algorithm α and t, where α controls the influence of the density of voxels on the manifold and t controls the scale of eigenvalues of the diffusion operator. According to previous recommendations, we used α at 0.5 and t at 0, preserving global relations between voxels in the embedding space (Margulies et al. 2016; Vos de Wael et al. 2020; Huang et al. 2023). Here, t = 0 indicated that the diffusion time was derived from an automated estimation using a damped regularization process (Margulies et al. 2016; Vos de Wael et al. 2020; Huang et al. 2023). To better interpret the spatial pattern of WMFC gradient, the gradient scores, which reflect the relative position of brain voxels along the gradient axes (Margulies et al. 2016), were projected into functional organizations.

In addition, to localize differences between BOLD and FDG WMFC gradients, a rank of each voxel within the WM mask was calculated, i.e. all voxels were labeled according to the gradient scores in ascending orders (Yang et al. 2021). Then, the difference was defined as that BOLD gradient rank minus FDG gradient rank.

Cell type-specific gene expressions association with WMFC gradients

The functional gradients are related to brain-wide gene expressions (Huntenburg et al. 2018; Sydnor et al. 2021). Thus, we next aimed to understand the underlying microscale contributions to the distinct and shared characteristics of the principal gradients between BOLD and FDG WMFC. The Allen Human Brain Atlas (AHBA) microarray expression data (Hawrylycz et al. 2012) bridged the gap between macroscale gradients of WMFC and microscale gene expression (Fornito et al. 2019). The AHBA microarray-based gene expression data provided high-resolution coverage of nearly the entire brain, with 3,702 spatially distinct tissue samples taken from 6 neurotypical postmortem brains. In line with our previous studies (Li et al. 2021c, 2021d), we first obtained the gene expression maps in WM from Neurosynth-Gene (https://www.neurosynth.org/genes/). Because spatially comprehensive maps of cell type densities across the WM were unavailable, we used previously defined cell-class gene sets (Seidlitz et al. 2020). We organized cell types into 7 canonical classes: microglia, endothelial cells, oligodendrocyte precursors, oligodendrocytes, astrocytes, and excitatory and inhibitory neurons (Seidlitz et al. 2020). The averaged gene expressions of cell type-related gene sets were then considered as the spatial patterns of cell types.

A multivariate linear regression model combining 7 cell types was used to explore the cellular contributions to the rank difference of principal gradients (G1) between BOLD and FDG WMFC, using the relaimpo (relative importance of regressors in linear models, version 2.2–5) in the R package. Relative importance metrics can be used to address linear regression with multiple collinear regressors (Groemping 2006). The model was defined as follows:

where |${\mathrm{Diff}}_i$| is the rank difference of the principal gradient score at a given voxel |$i$|⁠. After determining the percentage of explained variance in rank difference of the WMFC gradient from these predictors, the relative importance of each cell type regression contribution to the bootstrapping regressor model was assessed (Mandal et al. 2020).

Null models

The P-values in this study were tested against 2 categories of null models. The first null model was based on a blockwise permutation test for temporal correspondence. The second null model was based on Moran Spectral Randomization to minimize the effect of spatial autocorrelation for testing linear associations between BOLD and FDG WMFC topological properties (DC and gradients). The third null model was based on a permutation test for gradient rank difference. The detailed information is provided in Online-Only Data Supplement, Supplemental Methods.

Results

Temporal correspondence of the BOLD and FDG signals

We decomposed WM BOLD data into 10 spatial ICs to determine WM functional networks. Each IC represented voxels of a network that exhibited similar patterns of BOLD fluctuations over time. WM functional networks included: the auditory network (AN), the default mode network (DMN), the dorsal attention network (DAN), the visual association area (VAA), the sensorimotor network (SMN), the primary visual network (VN), the corpus callosum (CC), the brainstem network (BS), the sensorimotor-auditory network (S&A), and the cerebellum network (Cere; Fig. 1B).

We obtained MI values between BOLD and FDG signals for each participant in each WM functional network. We found that DMN (P = 0.047), VN (P = 0.001), and S&A (P = 0.003) exhibited significant MI values across all participants (Fig. 1b).

Local architecture of the BOLD and FDG WMFC

Little is known about the putative coupling between BOLD and FDG WMFC. DC is a reliable and common measure used to count the weights of connections on any voxel (Fig. 2A, bottom row). Based on visual inspection of group-averaged DCs, the 2 spatial distributions had some commonalities, where the splenium of CC and Cere exhibited high DC values. We then computed the Spearman's correlation coefficient between their DC values across different thresholds (r) from 0.0 to 0.6 (interval = 0.1). For example, BOLD DC spatially correlated with FDG DC at a threshold r = 0.6 (Rho = 0.81; P_Moran = 0.0002, FDR-corrected; Fig. 2B). We observed statistically significant correlations across all thresholds (all P_Moran < 0.05, FDR-corrected; Fig. 2C), suggesting a consistently similar local topology between BOLD and FDG WMFC. In addition, we validated the spatial similarities between the BOLD and DC at individual levels. We found significant positive correlations between 2 modalities at thresholds (r) ranging from 0.3 to 0.6 (Binomial test; all P_moran < 0.05, FDR-corrected; Fig. S2).

BOLD fluctuations concentrations with FDG WMFC

We computed the functional harmonics (i.e. the eigenvectors of the graph Laplacian) of the FDG WMFC (Fig. 3A). The functional harmonics exhibited a smoothly varying pattern across the WM between positive and negative polarities (Fig. 3B).

To quantify the coupling strength between the 2 functional modalities, we projected BOLD fluctuations on FDG functional harmonics; i.e. for each time point, the spatial pattern of BOLD activity was represented as a weighted linear combination of harmonic components of FDG WMFC. We divided the harmonics into 3 portions, i.e. low-, medium-, and high-frequency, with equally averaged energies. The aligned components with low frequency represented BOLD fluctuations that varied smoothly across the FDG WMFC at a single moment. The liberal components with high frequency denoted fluctuations that varied highly across the network (Medaglia et al. 2018).

The aligned components were concentrated within the AN, DMN, DAN, VN, CC, S&A, and Cere functional networks of the WM (all P_perm = 0.003, FDR-corrected; Fig. 3C and D). In addition, the liberal components were concentrated on the WM DMN (P_perm = 0.007, FDR-corrected), SMN (P_perm = 0.007, FDR-corrected), and BS network (P_perm = 0.007, FDR-corrected; Fig. 3E and F). Notably, we found that the WM DMN shared both aligned and liberal fluctuations, indicating that the content of BOLD in this WM system was complex to underlying glucose metabolic connectivity.

Global correspondence between BOLD and FDG WMFC gradients

We evaluated the BOLD and FDG WMFC gradients using the diffusion map-embedding technique (Fig. 4A), and we captured the principal gradient (G1) and the secondary gradient (G2) on the BOLD and FDG WMFC. The G1 and G2 were scattered in a triangular distribution pattern similar to a canonical distribution (Fig. 4B and C, middle panel). The BOLD WMFC G1 was anchored at 2 extreme ends involving the VN and SMN, and another involving the BS and Cere networks (Fig. 4B). The 2 ends of FDG WMFC G1 involved the VAA, SMN, BS, and Cere networks (Fig. 4C). The G1 pattern included the VN and AN on one end, and the DMN on another end, which was consistently distributed across the BOLD and FDG. We further determined that the G1 of BOLD WMFC spatially correlated with the G1 of FDG WMFC across whole WM voxels (Rho = 0.66, P_Moran = 0.008; Fig. 4D).

To determine distinct spatial variations, we directly compared the principal gradients (G1) between the BOLD and FDG WMFC using voxel ranks (Fig. 4E, left and middle panels). Voxel-wise rank comparisons revealed shifts in the VN and AN to the BOLD WMFC G1, and a converse shift in high-level association voxels to the FDG WMFC G1 (Fig. 4E). We found that the BOLD WMFC G1 dissociated from FDG WMFC G1 in the VN (P_perm = 0.0002, FDR-corrected) and AN (P_perm = 0.0002, FDR-corrected), whereas there was a converse dissociation of FDG WMFC G1 from BOLD WMFC G1 in the DMN (P_perm = 0.0002, FDR-corrected; Fig. 4F, right panel).

After obtaining gene expressions within each cell type, we used a multiple linear regression model combining all cellular factors (Fig. 5A). The model explained 13% of the variance in the rank difference of G1 between BOLD and FDG WMFC (F_{(7, 5613)} = 116.7; P < 2.2e−16), and glial cell type gene expressions significantly predicted the BOLD WMFC G1 (Fig. 5B; Table S1).

The reproducibility of correspondence between BOLD and FDG

Because hemodynamic response function (HRF) is essential in analyzing the BOLD-fMRI (Li et al. 2019), we re-evaluated the topological correspondence between BOLD and FDG WMFC. To minimize any HRF effects, we used a blind-deconvolution technique for BOLD (Wu et al. 2013, 2021). After reconstructing BOLD WMFC, we compared BOLD to FDG WMFC following the same parameters as previously mentioned. We also found a highly similar DC between the 2 modalities across a range of thresholds (all P_Moran < 0.05; Fig. S3). We then compared 10 ICs from deconvoluted BOLD data to BOLD data and found a similar spatial pattern of WM functional networks (Fig. S4). Finally, we re-evaluated the principal gradient of deconvoluted BOLD data. We found that the G1 of deconvoluted BOLD WMFC still correlated with the G1 of FDG WMFC (Rho = 0.67, P_Moran = 0.0008; Fig. S5).

Discussion

Using fPET-fMRI with a continuous infusion of FDG, we characterized the correspondence between intrinsic BOLD functional connectivity and glucose metabolic temporal connectivities within the human brain WM. Consistent with previous reports characterizing the direct relationship between BOLD signals and metabolic demand variations (Di et al. 2012; Wehrl et al. 2013; Amend et al. 2019; Jamadar et al. 2021; Guo et al. 2022; Voigt et al. 2022), we found that measurable functional networks from BOLD were accompanied by substantial correspondence to glucose metabolic temporal correlations. Specifically, we found that DMN shared aligned and liberal BOLD fluctuations based on FDG WMFC, indicating that the content of BOLD fluctuations in the WM DMN might involve a complex underlying FDG uptake. In addition, we found that functional gradients were dissociated between BOLD and FDG WMFC in the VN and DMN. Moreover, we found that cell type-specific gene expressions, especially of astrocyte cells, usually explained the rank differences of principal gradients between BOLD and FDG WMFC, indicating that glial cell gene expression patterns contributed to the difference between the 2 modalities. These findings may fill a gap in our understanding of the glucose metabolic processes of BOLD WMFC during the resting state.

Temporal synchronization of BOLD and FDG signals

Glucose is the primary energy substrate in the brain. [¹⁸F]FDG-PET provides the opportunity to characterize metabolic elements of brain connectivity based on glucose uptake (Yakushev et al. 2017; Jamadar et al. 2021). In contrast to BOLD-fMRI, [¹⁸F]FDG-PET quantifies cellular glucose consumption and can uncover the metabolic basis of neuronal and glial activities (Amend et al. 2019). The WM comprises more than 50% of the brain, including myelinated axons and glial cells. Compared with GM, metabolic demands in WM are mainly glial activities (Harris and Attwell 2012; Bonvento and Bolanos 2021; Xiang et al. 2021), such as oligodendrocyte cells, astrocytes, and microglia. In addition, as mentioned previously, several studies have successfully detected the functional information using fMRI (Gawryluk et al. 2014; Ji et al. 2017; Peer et al. 2017; Gore et al. 2019; Li et al. 2019, 2020a, 2020b; Guo et al. 2022). We were therefore concerned with the temporal synchronization between BOLD and FDG signals and the spatial correspondence between BOLD and metabolic connectivity.

For the temporal correspondence, we tried to link BOLD changes in WM to underlying dynamics of FDG glucose uptake using MI values. The MI has been shown to be robust in quantifying the relationship between any 2 waveforms, and it is not limited to specific function types. We found that BOLD signals shared more MI with FDG in WM DMN, VN, and S&A, indicating the associations of temporal correspondence between BOLD and FDG signals. These WM networks were identified by group-level sICA. One of the hallmarks of brain networks is their modular organization (Sporns and Betzel 2016). The modular organization is essential for both evolution and development, involving determining major building blocks or subnetworks that are internally densely connected and externally sparse (Clune et al. 2013). A previous study identified symmetrical WM functional networks, highly corresponding to GM networks at rest, highlighting their functional roles (Peer et al. 2017). Li et al. (2020c) have used ICA approach to the GM analysis of the fPET data set. However, we used BOLD data to obtain the WM functional networks consistent with previous studies in studying WM functional organizations (Peer et al. 2017; Basile et al. 2022).

BOLD and FDG WMFC architectures and their drivers

Recent simultaneous PET-fMRI studies have investigated the correspondence between glucose metabolism and BOLD functional networks in the human GM (Di et al. 2012; Savio et al. 2017; Jamadar et al. 2021; Voigt et al. 2022) and in rats (Wehrl et al. 2013; Amend et al. 2019; Ionescu et al. 2021), indicating the complementary strengths of FDG and BOLD in measuring the intrinsic connectivity of the brain. In contrast to these studies, this work investigated human brain WM function using brain connectivity analyses (Ji et al. 2017; Peer et al. 2017; Li et al. 2019, 2020a, 2020b, 2021d). Previous studies focused on WM functional networks derived from BOLD fluctuations driven by a still incompletely understood convolution of cerebral blood flow, cerebral blood volume, and cerebral metabolic rate of oxygen; herein, we studied the relationships between BOLD and FDG WMFC. Our results were similar to 2 recent simultaneous fPET-fMRI reports (Jamadar et al. 2021; Guo et al. 2022). One study aimed in part to relate WM BOLD fluctuations and connectivity to local metabolism, revealing that BOLD signals in WM were related to variations in metabolic demand (Guo et al. 2022). The other study applied the functional connectome method to highlight the complementary strengths of BOLD and FDG in measuring the intrinsic connectivity in GM (Jamadar et al. 2021). However, in contrast to these studies (Jamadar et al. 2021; Guo et al. 2022), our focus was on characterizing the shared and different multiscale topological properties of BOLD and FDG WMFC. DC is an example of a strictly local measure; it characterizes only a single node or a single voxel (Takeuchi et al. 2015). It is important to note that while BOLD and FDG data record different neurophysiological levels, they both record from the brain WM, assaying the same underlying WMFC networks, and their WMFCs show strong relationships with DC across a range of thresholds, suggesting an underlying synchrony between hemodynamic processes and FDG uptake. In addition, we validated the correspondence between BOLD and FDG WMFC after considering the HRF impact on BOLD data, revealing the reproducibility of underlying synchrony between hemodynamic processes and FDG uptake.

The BOLD activity underlying FDG WMFC

We added a critical perspective for broader interests in the constraint of FDG uptake to hemodynamic responses. By definition, any pattern of brain activity can be expressed using this functional basis as a superposition of functional harmonics (Glomb et al. 2021). We thus used the functional harmonics of FDG WMFC to decompose BOLD fluctuations. Using the harmonic decomposition of FDG WMFC, we first showed that observed brain activity was preferentially expressed using low-frequency components (graphs). These activity patterns fit better to the connectome constraints, indicating that activity patterns expressed smoothness on the connectome (Preti and Van De Ville 2019). The BOLD signal was filtered into 2 parts based on functional harmonics of FDG WMFC: one was aligned components; i.e. keeping low-frequency components coupled with the FDG WMFC; the other was liberal components; i.e. keeping high-frequency components decoupled with the FDG WMFC. As a complement to previous findings implicating spatial correspondence of BOLD and FDG architectures, our results indicated the importance of FDG WMFC organization, and the central role of the DMN in this structure. The content of WM BOLD fluctuations in the DMN for underlying metabolic connectivities was complex. Aligned components within the DMN implicated coupling between the 2 modalities. In contrast, dissociations between BOLD fluctuations and FDG uptake within the WM revealed functional heterogeneity in this network (Stiernman et al. 2021). A question is why do some functional networks appear to have liberal components? One possibility is that both the BOLD and FDG systems carry some unique neurophysiological information. Compared with BOLD, FDG represents a more direct index for measuring neural activity independent of neurovascular coupling (Jamadar et al. 2020). More broadly, the liberal components re-emphasize the complex organization of the WM, in which voxels are linked to others by different connectivity modalities.

Topological gradients of BOLD and FDG WMFC

Brain fundamental connection features pertain to different brain modalities, such as structure, correlated activity, or gene expression. Recent advances in mapping brain areas provide the basis for investigating the significance of their spatial arrangements, which are considered a core organizational component of human brain structures (Huntenburg et al. 2018). In the present study, we described the dominant gradients in BOLD, and FDG WMFC features that spanned from the VAA and sensorimotor regions to the Cere network. The Cere has traditionally been considered to have extensive connectivity with motor aspects of extracerebellar structures. In contrast, recent neuroimaging studies have reported that the human Cere was involved in motor control and in cognitive control processing. The central axis of cerebellar functional regions follows a gradual organization that progresses from primary (motor) to transmodal (DMN) regions (Guell et al. 2018). The principal axis of BOLD and FDG WMFC is like a meta-networking information of complex functions, i.e. from low-level to high-level networks, and then to transient meta-networks (Herbet and Duffau 2020), involving highly local networks, which process visual, sensory, and auditory inputs. Widely spread functional networks sustain high-level cognitive and affective processing in the brain. The brain integrates information from specialized networks for more complex and flexible cognition. The similarity of the principal gradients of the 2 modalities provides an understanding of how the spectrum of WM function relates to brain glucose metabolism.

Notably, we also observed specific hierarchical organization between BOLD and FDG WMFC gradients. The BOLD WMFC gradient exhibited a dissociation from FDG uptake in the VN. In contrast, FDG WMFC gradients displayed a converse dissociation from the hemodynamic response in the DMN, reflecting the segregation of anterior–posterior communities between the 2 modalities. This result is consistent with previous findings, which reported that metabolism in widespread DMN remained high during both rest and tasks, despite negative BOLD responses, indicating functional heterogeneities in this network (Stiernman et al. 2021). In addition, we found a complex content of BOLD in the DMN regarding the underlying FDG uptake. In contrast, the DMN exhibited dissociation between the 2 modalities from global topological properties, indicating the distinct information integration in this network in different context-sensitive manners, i.e. as a function of current cognitive demands (Herbet and Duffau 2020). Together, the DMN was identified as a core network in both hemodynamic responses and glucose metabolism.

Simple models of complex WM functional networks

Understanding the rules that guide the principal gradients of WMFC is a critical cognitive neuroscience goal for WM functional studies. We identified a set of cell types to explain functional variation differences between the BOLD and FDG WMFC across WM voxels. Our findings suggest that spatial gene expressions of glial cells had a relatively high predictive capacity in globally determining WM functional organizations differences between hemodynamic responses and FDG uptake/metabolism, indicating that glial cells may influence information processing for different modalities. Functional gradient shifts have been reported to occur throughout development (Dong et al. 2021), and early human brain developmental studies from single-cell analysis highlight the heterogeneity of human cells (Eze et al. 2021). Collectively, these investigations build on past findings aimed to identify the role(s) of WM functions derived from BOLD, and furthermore reveal a similar microscopic architecture of WM functional organizations from different modalities.

Methodological considerations

To improve our understanding of the associations between neurohemodynamics and metabolism in WM, future considerations may be addressed. One particular concern is that slow changes in systemic physiology may lead to a portion of WM functional signals. Thus, future studies investigating BOLD signals should account for multiple physiological measures (Ozbay et al. 2018), such as an MRI scanning system measuring heart rate, respiratory rate and depth, and photoplethysmography or near-infrared spectroscopy recording peripheral measurements of vascular tone. Second, we tried to link microscale gene expression of cell types to the macroscale gradient differences between BOLD and FDG WMFC. We found that gene expressions of glial cells were associated with spatial distribution. However, gene expression data was available from AHBA data set, probably limiting the explanation of the contributions. Third, the samples in this study were comprised of young healthy adults (ages ranging from 18–23 years), limiting the understanding of glucose metabolic mechanisms of BOLD WMFC for other populations. Future studies should analyze a broader population base to investigate the relationships between BOLD and FDG WMFC patterns. Finally, we have presented the temporal and spatial correspondence between BOLD and FDG signals, trying to illustrate the shared and specific functional information to decode brain function and brain WM networks. However, the shared information may be caused by the shared anatomical information in the 2 modalities, such as heterogenous cell and vessel densities across WM. Future studies should investigate the effects of anatomical information on BOLD and FDG images to better understand the correspondence between the 2 measures.

Conclusion

In summary, we described the correspondence between hemodynamics and underlying metabolic dynamics within the WM. BOLD and FDG data display shared and complementary views of the WMFC from the spatiotemporal topological architectures. Notably, cell type-specific gene expressions may be responsible for the principal gradient differences between BOLD and FDG WMFC. Collectively, our findings add to the current understanding of WM BOLD functions and provide insights into the coordination of hemodynamics and energy metabolism.

Author contributions

W.L. and H.C. led the project. J.L., W.L. and H.C. were responsible for the study concept and the design of the study. J.L. and W.L. created the figures and wrote the manuscript. G-R.W. performed PET data analysis. M.S. and J.X. analyzed the neuroimaging data. Y.M. and S.Y. performed gradient analysis. All authors reviewed and commented on the manuscript.

Acknowledgments

The authors are grateful to all participants in this study and thank Jamadar et al. published the simultaneous [¹⁸F]FDG-fPET and BOLD-fMRI data in OpenNeuro.

CRediT author statement

Jiao Li (Conceptualization, Investigation, Writing—original draft, Writing—review & editing), Guo-Rong Wu (Methodology, Writing—review & editing), Mengyuan Shi (Methodology, Writing—review & editing), Jie Xia (Methodology, Writing—review & editing), Yao Meng (Validation, Writing—review & editing), Siqi Yang (Validation, Writing—review & editing), Huafu Chen (Funding acquisition, Writing—review & editing), and Wei Liao (Conceptualization, Writing—original draft, Writing—review & editing)

Funding

The National Science and Technology Innovation 2030 Major Program (2021ZD0201701); National Key Project of Research and Development of Ministry of Science and Technology (2022YFC2009900 and 2022YFC2009906); National Natural Science Foundation of China (82202250, 62276051, 62036003, 61876156); China Postdoctoral Science Foundation (BX2021057, 2022 M710615).

Conflict of interest statement: None declared.

Data availability

Simultaneous [¹⁸F]FDG-fPET and BOLD-fMRI data were openly available on OpenNeuro (accession number: ds002898) provided by Jamadar et al. (2020).

References