The complexity of the relationship between spontaneous brain activity and glucose metabolism

Brain glucose metabolism as assessed by [ 18 F]FDG positron emission tomography (PET) is expected 37 to be significantly related to resting-state functional MRI (rs-fMRI) activity and functional 38 connectivity (FC), but the underlying coupling model is still incompletely understood. Employing 39 simultaneous acquisitions, we related [ 18 F]FDG standard uptake value ratio (SUVR) to 50 features 40 pertaining to rs-fMRI 1) signal, 2) hemodynamic response, 3) static and 4) time-varying FC, and 5) 41 phase synchronization. To assess which rs-fMRI variables better describe SUVR across regions, we 42 employed a hierarchical approach, identifying the model at population level, and then estimating it 43 on individual data. Multilevel modelling explained around 40% of the SUVR variance, with signal- 44 related features as the most relevant fMRI variables. When the model was used to characterize 45 between-network variability of the SUVR-fMRI coupling, the ranking changed. We demonstrate that 46 local activity and synchronization are the most important predictors of glucose metabolism, while 47 large-scale FC properties gain importance within specific networks.


Introduction
1. which is the strength of the bivariate association between these rs-fMRI features and SUVR across 123 the whole brain? And then, since regions with high vs. low metabolic consumption are expected to 124 have quite different structural and functional properties 1,3,4 , does this coupling change according to 125 the ranking of brain nodes based on SUVR? 126 2. is it possible to explain group level SUVR variance across regions by combining rs-fMRI features, 127 for the first time, into a multiple regression model? Is the group of selected features more populated 128 by local or large-scale brain network metrics, and does it account for between-subject variability 129 (BSV) 30 ? Finally, which of the previously identified rs-fMRI features are more important to explain 130 SUVR when multilevel modelling is performed across RSNs, i.e., which is the between-network 131 variability (BNV) of the SUVR-fMRI association? 132 133

Feature extraction and preliminary evaluation of rs-fMRI variables 143
The flowchart describing the preprocessing and preliminary analysis of the [ 18 F]FDG PET and rs-144 fMRI data is shown in Figure 1 (see the Methods section for details). 145 The [ 18 F]FDG PET variable of interest is the SUVR, which was extracted for every region of the 146 Schaefer cortical atlas 32 (200 parcels, supplemented by 18 subcortical regions 33 ) in each subject, and 147 will be considered as the dependent variable in every modelling approach from here onward. 148 The 50 rs-fMRI variables, extracted at the single-subject level and a priori subdivided into 5 pools, 149 are reported in Table 1: the signal pool (1) contains features related to the basic statistics of the BOLD 150 time series (median, variance, skewness), its complexity, its low-frequency fluctuations (ALFF), local 151 coherence (ReHo) and high-amplitude events (peaks-BOLD); in the HRF pool (2), then, we placed 152 the amplitude of the HRF peak, calculated using a blind deconvolution method 34 , and the HRF 153 correlation structure across regions described by means of graph properties; the sFC pool (3) 154 characterizes FC calculated across the entire fMRI scan with graph theory metrics; the tvFC pool (4) 155 assesses graph metrics' temporal variability across sliding windows 13 ; finally, the PC pool (5) 156 characterizes FC as coherence of BOLD phase 29 . Both rs-fMRI time series and [ 18 F]FDG SUVR data were parceled using the Schaefer cortical atlas 159 (200 ROIs) and 18 subcortical ROIs. The parcel-wise rs-fMRI data were used to extract fifty features 160 representative of five pools, i.e., 1) signal, 2) HRF, 3) sFC, 4) tvFC, 5) PC. The PET-fMRI coupling 161 was investigated using bivariate correlation and multivariable multilevel modelling across subjects 162 and across fMRI-based RSNs. 163 The Spearman's correlation matrix between the 50 rs-fMRI variables at group median level (i.e., by 168 taking the parcel-wise median value of each feature across subjects) was computed (Figure 2a), in 169 order to assess the relationships between the extracted features and their degree of redundancy: the 170 clustering into 5 pools provided by a priori knowledge was fairly consistent with the observed 171 correlation structure, with signal, HRF and sFC features (upper block) being clearly distinguished 172 from tvFC features (lower block), which they are negatively correlated with, and PC variables 173 demonstrating lower correlation with the rest. However, it was also noticeable that strong correlations 174 between many variables were present, especially for the tvFC pool, and that a feature selection step 175 was going to be necessary to use these variables in a numerically sound multivariable model of 176 SUVR: the condition number ( ), which quantifies the level of correlation between predictors in a 177 multiple regression context (i.e., their multicollinearity), was high ( ( ) = 70.58), way beyond the 178 acceptability range 35 , and this is known to result in unstable and unreliable models (see Methods). 179 180 181

SUVR vs. rs-fMRI: bivariate relationships 182
Before moving to the multiple regression framework, we began by investigating bivariate associations 183 between SUVR and the extracted rs-fMRI variables at the group level, in the so-called naïve average 184 data approach (NAD), as done by many previous studies 23,25,26 ; here, however, a much wider range 185 of fMRI-derived variables was explored. Many significant spatial associations between SUVR and 186 rs-fMRI features were detected across the 218 analyzed regions, as assessed through Spearman's rank 187 correlation (p = 0.05 significance level) with false discovery rate (FDR) multiple comparison 188 correction 36 . The correlation coefficients are reported in Figure 2b.  The pattern of Spearman's correlations (FDR-corrected, non-significant values shown in  193  white) among rs-fMRI features, assessed at the group level and divided according to the pool to which  194 they have been assigned (1) signal, 2) HRF, 3) sFC, 4) tvFC, 5) PC), is shown in (a). The rs-fMRI  195  features are tested for association with group median SUVR across 218 brain regions (b) via  196 Spearman's correlations (significant values after FDR correction are indicated with an asterisk). 197

199
The strongest positive associations were with 1) ReHo ( = 0.45, p < 0.001), 2) s-BC ( = 0.4, p < 200 0.001), and 3) SampEn-BC ( = 0.44, p < 0.001), i.e., respectively 1) a measure of local 201 synchronization of BOLD, 2) a sFC graph metric, betweenness centrality (BC), which describes a 202 node in terms of its global connections in a graph, and 3) a measure of temporal complexity of the 203 BC time series. The strongest negative correlations were mdiff-BC ( = -0.42, p < 0.001) and CV-204 BC ( = -0.42, p < 0.001) in the tvFC pool, both measures of temporal variability of BC. 205 In general, it can be noted that positive associations emerged for the majority of the signal-based, 206 HRF and sFC-related features, while tvFC metrics displayed a consistent and never previously 207 reported negative association with SUVR (Figure 2b) Overall, it can be noted that the detected spatial correlations were at best moderate. 215 216 SUVR-fMRI associations are strengthened in low SUVR nodes 217 As the relationship between [ 18 F]FDG PET and rs-fMRI could be spatially heterogeneous, 218 Spearman's correlations were also re-evaluated across groups of nodes selected according to linearly 219 increasing percentiles of the SUVR distribution, i.e., by progressively retaining the highest SUVR 220 values, from the 1 st up to the 85 th percentile, as well as linearly decreasing percentiles of the same 221 distribution, i.e., by retaining the lowest SUVR values, from the 100 th to the 15 th (Supplementary 222 Figure S1). The purpose of the analysis was to verify whether SUVR-fMRI associations would be 223 strengthened in high SUVR nodes or, conversely, in low SUVR nodes, as SUVR provides a clear 224 ranking of brain regions, that is expected to be related to crucial structural and functional properties, 225 e.g., neuron-to-glia ratio, richness in neuroreceptors, excitatory-inhibitory activity 1,3,4 . 226  We then identified the percentile threshold corresponding to the highest total correlation value across 244 features: the spatial pattern of the 87 nodes below the 40 th percentile of the SUVR distribution is 245 shown in Figure 3b. These parcels, where the SUVR-fMRI association is emphasized, mainly belong 246 to temporal/limbic areas (including hippocampus), sensorimotor cortices, and subcortical regions, 247 such as cerebellum and globus pallidus (Supplementary Figure S1). We then set out to assess which combination of rs-fMRI features was better able to explain SUVR 256 across brain regions with multiple regression and multilevel modelling, in a fully data-driven way. In 257 multilevel modelling, the model structure is usually known, or selected at the lower level, i.e., at the 258 individual level 30 . However, as significant BSV in the SUVR-fMRI association is expected, 259 especially for the rs-fMRI features, we chose to identify the model predictors at the population level 260 (NAD approach), thus exploiting the denoising properties of averaging. The model structure selected 261 at the group median level was then used for multilevel modelling across subjects to characterize the 262 BSV of the SUVR-fMRI association, as shown in (Figure 4, top), trying to fully capitalize on the 263 fact that [ 18 F]FDG and rs-fMRI data were acquired in the same subjects.   In this case we found that the BSV in the SUVR-fMRI association is clearly non-negligible. In 320 fact, the group-level ! estimates are very close to those obtained using the NAD approach, 321 confirming the adequacy of the average approach in describing the relationship between the variables. 322 However, as expected, the R 2 of the overall model, i.e., considering BSV, was lower and equal to 323 0.245, due to the capability of the multilevel mixed-effect approach to keep into account both between-and within-subject variability. The R 2 values of the subject-level models are reported in 325 Supplementary Figure S7a, and they display high variability (from 0.05 to 0.45). 326 The median across subjects of the model's residuals !" , which highlight how well the SUVR of each 327 region is explained by the identified model, can be visualized in Figure 5c. Notably, high positive 328 values are present in posteromedial cortex (posterior cingulate cortex (PCC) in particular) and 329 subcortex (putamen): these areas identify nodes with high SUVR values which are not satisfactorily 330 explained by the available rs-fMRI features. Importantly, this deficiency in explanatory power is 331 highly consistent across subjects, as evidenced by the low variability (CVs%) of the residuals in those 332 areas (Figure 5d). 333 errors for the fixed effects ! , which represent the parameters that best explain SUVR across regions 337 at the group level (a). The relative importance weights produced by dominance analysis (DA), 338 highlighting the proportion of the multivariable multilevel model R 2 explained by each predictor 339 (general dominance) (b). Across-subject median (c) and CVs% (d) of weighted residuals !" of the 340 multilevel model, plotted on the brain cortex and subcortex. 341 342 343 SUVR vs. rs-fMRI: multilevel model across networks 345 Finally, the log-linear model with the 9 selected rs-fMRI predictors was tested in its ability to describe 346 the expected network-level variability of the SUVR-fMRI association (BNV), i.e., by considering 347 only the parcels referring to specific brain networks. RSNs were here grouped according to the 348 Schaefer' functional atlas in its 17-RSN partition 32 , supplemented by 18 subcortical anatomical 349 regions 33 . A multilevel modelling framework was again employed, but with RSNs as the grouping 350 factor for individual level data, instead of subjects (Figure 4, bottom). 351 The fixed effects # and their SEs for the between-network model are reported in (Figure 6a) (Figure 6b). 359 Notably, the R 2 of model prediction considering network-wise estimates is markedly lower than when 360 subjects are used to cluster nodes. As shown in (Figure 6c), the single RSNs are highly heterogeneous 361 in terms of model R 2 , ranging from around 0 to 0.32, with an overall prediction with R 2 = 0.147. The BNV of the SUVR-fMRI association is measured by the random effects #$ for each network, 373 with some RSNs displaying significant distance from the model estimates # of the "average-374 network". To better assess this variability, the nine rs-fMRI predictors' parameter estimates #$ (i.e., 375 sum of fixed effects # and random effects #$ for every network j) were plotted (Figure 7a). We 376 can observe that most predictors included in the multivariable model display heterogeneity across 377 networks in their relationship with SUVR, with either positive or negative associations depending on 378 the specific RSN, which cannot be captured by the average situation described by the fixed effects 379 # of (Figure 6a). Some predictors show notably consistent spatial patterns, and therefore, to assess 380 their similarity, the correlation between their random effects was evaluated across networks ( Figure  381   Individual network parameter estimates ( #$ , sum of fixed effects # and random effects #$ , which 397 describe the variability from the fixed effect for each RSN j), plotted on the brain surface for each 398 predictor (a). Correlation matrix (non-significant values in white, p = 0.05, uncorrected) of the nine 399 predictors' random effects #$ across RSNs (b). Correlation matrix (non-significant values in white, 400 p = 0.05, uncorrected) of the RSNs' parameter estimates ( #$ = # + #$ ) across predictors (c). 401 402

404
In this work, we thoroughly investigated the spatial relationship between a wide range of features 405 extracted from rs-fMRI and simultaneously acquired [ 18 F]FDG PET, while also accounting for the 406 variability across subjects (i.e., BSV) and networks (i.e., BNV) in this relationship. 407 408

Relationship between SUVR and rs-fMRI through bivariate correlation 409
In addition to the rs-fMRI variables that have already been associated to SUVR, i.e., ALFF, ReHo, 410 sFC strength 22,23,25 , we extended our assessment to a wide variety of previously unexplored features, 411 such as time-varying functional connectivity (tvFC) and HRF-related variables. 412 To our knowledge, the relationship between [ 18 F]FDG metabolism and FC temporal variability has 413 never been tested before. While it is established that regions with stronger static FC tend to have 414 higher cerebral blood flow 41 and higher energy metabolism 25 , possibly reflecting the fact that they are 415 also more strongly connected anatomically 42 (with recent work highlighting that structural 416 connectivity graph properties are positively associated with [ 18 F]FDG SUVR 43,44 ), the tvFC coupling 417 with glucose metabolism is not established. 418 We found that tvFC (as captured by graph theory metrics' temporal variability) has a moderate-strong 419 negative association with SUVR. The interpretation of this finding can be supported by knowing that 420 sFC and tvFC graph metrics are negatively correlated, as clearly shown by the correlation matrix 421 among rs-fMRI predictors (Figure 2a)  Additionally, a strong negative relationship is found between SUVR and the number of BOLD 433 pseudo-events (peaks-BOLD), which is related to the interpretation of the BOLD signal as a point 434 process, with sparse neural events governing its dynamics 48 . While puzzling at first, one interpretation 435 might come from considering that higher local oxygen consumption by active neurons is associated 436 with decreased positive BOLD fluctuations 49 , and therefore the higher the number of BOLD peaks 437 and extreme events, the lower the oxidative metabolism and SUVR might be in that region.
We then examined how these relationships would be modulated by selecting parcels according to 439 their ranking in the SUVR distribution. In order to better probe the spatial relationship between SUVR 440 and rs-fMRI, which is heterogeneous across the brain 23,28 , we chose to explore the changes in 441 correlations selecting nodes from the SUVR standpoint, instead of according to FC properties 24,25 . 442 Interestingly, nodes with progressively higher SUVR, which are expected to be the richest in terms 443 of receptor density, local activity and inter-regional communications 1,3,4 , did not show different 444 relationships between metabolism and rs-fMRI, but the correlations with most rs-fMRI features 445 became significantly stronger when considering nodes with progressively lower SUVR (Figure 3). To our knowledge, this is the first study to investigate the [ 18 F]FDG-fMRI coupling using a 461 multivariable approach, attempting to identify the best subset of metrics, among a wide range of 462 fMRI-derived variables, to explain SUVR variability across regions. Moreover, to fully capitalize on 463 the fact that PET and fMRI data were acquired in the same subjects, we employed a multilevel 464 modelling approach, with the selection of the best features performed at the group (higher) level, and 465 modelling performed at the individual (lower) level, to characterize the between-subject variability 466 of the SUVR-fMRI association (Figure 4, top). The selected model consisted of nine rs-fMRI 467 variables (Figure 5 a, b) which represented all pools of features: signal (ApEn-BOLD, rApEn-468 BOLD, peaks-BOLD, ReHo, CV-ReHo), HRF (hrf-LE), sFC (s-BC), tvFC (CV-BC), PC (med-Leig). 469 The strongest predictors are related to the BOLD signal and its local synchronization properties 470 (peaks-BOLD, ReHo), which consistently emerged as relevant across all feature selection methods 471 (Supplementary Results). The fact that the SUVR-fMRI spatial coupling is emphasized when local 472 BOLD variables are involved might reflect the interplay between excitatory and inhibitory neural 473 populations 51 , which regulate CBF, a main ingredient in many fMRI-related features 17,18,52,53 . 474 Overall, the explanatory power provided by BOLD rs-fMRI reached a 40% of the SUVR variance at 475 the group level (24% across subjects). Zones of polarization in the model residuals emerged in 476 subcortical (putamen), posteromedial (PCC), and lateral frontal regions, which could mainly be 477 attributed to outliers with higher metabolism (Figure 5c), which are poorly explained by the available 478 rs-fMRI features in a consistent manner across subjects (Figure 5d). These results point to the idea 479 that the BOLD signal and FC, even though related to CBV, CMRO2 and CBF 17,18 , reflect the 480 metabolic architecture established by [ 18 F]FDG SUVR only partially, even in simultaneously 481 acquisitions, and that rs-fMRI FC and its graph metrics cannot be considered a proxy of glucose 482 metabolism. 483 Moreover, the individual model R 2 values were variable across subjects, highlighting the fact that the 484 SUVR-fMRI relationship displays significant between-subject variability, with subjects whose 485 BOLD signal and FC architecture are more related to SUVR, and others that have hardly any 486 relationship. 487 Next, we used the multilevel modelling approach and the identified predictors to characterize 488 between-network variability, exploiting the fMRI-derived RSNs to group the individual data in a 489 network-by-network fashion (Figure 4, bottom). 490 Importantly, the rs-fMRI predictors selected for the between-subject model proved to still be relevant 491 for evaluating the between-network SUVR-fMRI association, but their ranking, as assessed by 492 dominance analysis, changed noticeably (Figure 6b), with static and dynamic FC features (CV-BC 493 in particular) gaining importance in the model. 494 Moreover, when the network-wise effects are considered, significant positive and negative 495 associations emerge for each of the nine predictors (Figure 7a). These patterns of predictors have 496 some degree of similarity across networks, with a cluster of RSNs (subcortical, limbic, salience, 497 dorsal attention etc.) being highly correlated, which implies they have a similar SUVR-fMRI 498 multivariable association pattern (Figure 7c). Other networks, instead, seem to be more isolated in 499 their SUVR-fMRI coupling (default mode, visual, somatomotor etc.), possibly reflecting their 500 enrichment in high SUVR nodes (DMN, VIS) that are more difficult to explain using fMRI features 501 (Figure 5c). 502 These findings add to and enrich previous work highlighting between-network variability in the 503 SUVR-fMRI association through bivariate associations 23,28 . 504 Notably, the regions where SUVR-fMRI bivariate correlations are higher (Figure 3) seem to fall into 505 networks with fairly high R 2 values in the multivariable model (Figure 6c), confirming that the 506 SUVR-fMRI coupling is emphasized in these regions. It is also important to underline that marked 507 regional heterogeneity has also been described in the coupling between [ 18 F]FDG SUVR and local 508 CBF 54 , which is likely to underlie some of the variability detected here due to CBF contributions to 509 i.e., K1 (ml/cm 3 /min), describing tracer uptake through the blood-brain barrier, k2 (min -1 ), describing 526 its efflux into the venous blood, and k3 (min -1 ), quantifying the phosphorylation rate of the hexokinase 527 in neurons and glia 1,2 . There is therefore the possibility that some of the contribution of CBF to the 528 [ 18 F]FDG-fMRI coupling comes from SUVR, which is in fact highly correlated with the early, CBF-529 related frames of [ 11 C]PiB PET in healthy controls 59 , as well as with PET-derived CBF estimates 54 . 530 It is therefore likely that PET kinetic modelling will help disentangle the biological processes 531 underlying both BOLD rs-fMRI and static PET estimates. 532 533

Conclusion 534
In conclusion, we thoroughly investigated for the first time the spatial relationship between [ 18 F]FDG 535 SUVR and a wide range of features derived from rs-fMRI (pooled into 1) signal, 2) HRF, 3) sFC, 4) 536 tvFC and 5) PC-based features) using simultaneous PET/fMRI data. Selection of low SUVR parcels 537 led to a strengthening of SUVR-fMRI associations, implying the presence of a nonlinear relationship 538 for many features. Moreover, a novel multivariable multilevel modelling framework was employed 539 to identify the best subset of rs-fMRI predictors able to explain regional SUVR variance, highlighting 540 that predictors based on the BOLD signal and its local synchronization (ReHo and BOLD pseudo-541 events, in particular) are the ones that are more tightly related to [ 18 F]FDG SUVR across brain regions. This suggests a local contribution of CBF that should be tested for and, possibly, regressed 543 out from the BOLD signal. 544 Notably, the overall explanatory power provided by rs-fMRI on the regional metabolic variability did 545 not exceed 40% of the variance at the group level, with significant variability across subjects. When 546 multilevel modelling of the SUVR-fMRI coupling was carried out across networks, the selected 547 predictors were still relevant for description of RSN metabolism, but noticeable variability across 548 networks was present: new positive and negative associations emerged, and sFC and tvFC network 549 features gained importance. In conclusion, SUVR variability across parcels is only partly expression 550 of brain network organization described by rs-fMRI. 551 552 553 The dataset includes 26 healthy subjects from two studies: 11 subjects (8 males; 52.2 ± 10.4 years), 605 hereby referred to as dataset A (Munich) 23 , and 15 subjects (6 males; 64.7 ± 7.9 years), i.e., dataset B 606 All subjects were identically pre-processed to obtain local metabolism information from [ 18 F]FDG 631 PET data, and BOLD-based measures from rs-fMRI data, employing a pipeline similar to the Human 632 Connectome Project (HCP) minimal preprocessing pipeline 60 with the addition of PET processing. 633