To systematically analyze the spatiotemporal relationship between resting-state electrophysiology and fMRI signals, we simultaneously recorded the whole-brain rsfMRI and electrophysiology signals in the M1 and ACC in rats in both light-sedation and awake states (Fig. 1A). The placement of implanted electrodes in the M1 and ACC was verified using T2-weighted structural images (Fig. S1). Raw electrophysiology data were initially preprocessed to remove MR artifacts using a template regression approach (Tu and Zhang, 2022). An example of denoised LFP is shown in Fig. 1A and Fig. S2.
LFP and rsfMRI signals derive highly consistent RSN spatial patterns in lightly sedated rats
We first examined the spatial relationship between brain-wide rsfMRI signals and frequency band-specific LFP powers in lightly sedated rats. For each recording site, its associated RSN is obtained as the seedmap, calculated by voxel-wise correlating the regionally averaged BOLD time series of the seed (M1 in Fig. 1F and ACC in Fig. 2B) with BOLD time series of individual brain voxels. This map is conventionally defined as the resting-state functional connectivity (RSFC) pattern for the seed. To determine to what extent this RSN can also be derived using the LFP signal recorded from the same location (M1 or ACC), the band-limited LFP power was generated using a conventional LFP band definition (delta: 1–4 Hz, theta: 4–7 Hz, alpha: 7–13 Hz, beta: 13–30 Hz, gamma: 40–100 Hz). Figures 1B&C show cross correlations of the LFP power and BOLD signal in the M1 across the LFP spectrum (Fig. 1B, 1-Hz band interval) and for individual LFP bands (Fig. 1C), indicating that gamma band power is positively correlated with the BOLD signal, while lower-frequency bands display negative peak correlations with the BOLD signal. The lag of the BOLD signal is approximately 2 sec for all bands, consistent with the hemodynamic response function (HRF) previously reported in rodents (Scholvinck et al., 2010; Winder et al., 2017). Therefore, for each LFP band, we convolved its power with an rodent-specific HRF (Fig. 1D, (Tong et al., 2019)) to generate the LFP band-predicted BOLD signal, which was then voxel-wise correlated with the brain-wide rsfMRI signal, yielding the LFP band-derived correlation map (Fig. 1E, Fig. 2A).
We found that the gamma band power-derived correlation map is spatially highly consistent with the corresponding seedmap. For M1, the voxel-to-voxel Pearson correlation coefficient (CC) between the mean gamma-derived correlation map (Fig. 1E) and mean M1 seedmap (Fig. 1F) is 0.95 (Fig. 1K, R2 = 0.90), indicating 90% of the variance in the M1 seedmap can be explained by the gamma-derived map. In contrast, the spatial maps generated by lower-frequency bands are inversely correlated with the M1 seedmap, with a trend that CC being more negative as the frequency range of the band is lower (Figs. 1E-J, delta: CC = -0.78; theta: CC = -0.78; alpha: CC = -0.5; beta: CC = -0.34), consistent with the LFP-BOLD cross correlations for these bands shown in Fig. 1B-C. These relationships are repeatable in the ACC (Figs. 2A-G, delta: CC = -0.62; theta: CC = -0.3; alpha: CC = -0.12; beta: CC = 0.12; gamma: CC = 0.85). These data indicate that the spatial patterns of BOLD-based RSNs can be reliably obtained using the band-specific LFP signal.
To examine whether these findings happen to result from the specific frequency cutoffs we adopted for any LFP band, we repeated the analysis for all individual 1-Hz bands in the full LFP spectrum. Our data again show a gradual transition from negative to positive spatial correlations in LFP-derived correlation maps with the corresponding seedmap as the LFP signal increases from low to high frequencies (Fig. 1L for M1; Fig. 2H for ACC). Furthermore, as a control analysis, the gamma band power in the ACC was temporally shuffled, convolved with the HRF and the correlation map was recalculated (Fig. S3). In this case the spatial pattern observed in Fig. 2A disappeared, suggesting the LFP-derived spatial patterns are indeed specifically related to the LFP signal, instead of an artifact of HRF. We also confirmed that all our results are not sensitive to the rsfMRI data preprocessing step of global signal regression (Figs. S4-S5).
These data collectively indicate the same RSNs can be reliably derived using LFP power in lightly sedated rats, suggesting the critical involvement of neural activity in RSN spatial patterns.
Temporal correlations between LFP power and local rsfMRI signal are significant but considerably weaker
Given the high reliability of the gamma power in determining the spatial patterns of BOLD-based RSNs, it is tempting to assume the HRF-convolved gamma power ought to be able to reliably predict the rsfMRI time series from the same location (Fig. 3A&C). Surprisingly, our data show that the local rsfMRI signal displays considerably weaker temporal correlations with LFP powers. In the M1, the LFP-BOLD temporal correlations gradually change from negative to positive as the LFP signal goes from low to high frequencies—the same trend observed in spatial correlations (Figs. 1E-K), with, however, significant but considerably lower absolute magnitude in CC values (one sample t-tests, delta: CC = -0.20, p = 2.1 × 10− 57; theta: CC = -0.19, p = 7.1 × 10− 62; alpha: CC = -0.11, p = 1.2 × 10− 42; beta: CC = -0.06, p = 8.4 × 10− 15; gamma: CC = 0.37; p = 2.1 × 10− 58; number of scans = 159). The same results are observed in the ACC (one sample t-tests; delta: CC = -0.13, p = 1.9 × 10− 28; theta: CC = -0.04, p = 8.4 × 10− 7; alpha: CC = -0.03, p = 6.0 × 10− 5; beta: CC = -0.01, p = 0.1; gamma: CC = 0.18, p = 6.7 × 10− 42; number of scans = 172).
Notably, one difference in our calculation of spatial and temporal correlations may contribute to the disparity in their CC values. In our calculation of spatial correlations, we first scan-wise generated LFP- and BOLD-derived RSN maps, and these maps were then averaged within each group before the spatial correlations were calculated using two averaged maps (Figs. 1–2). However, the temporal correlations were first calculated for individual scans (Figs. 3A&C) before the resulting correlation values were averaged across scans, as averaging time courses first would diminish the actual signal due to the semi-random nature of spontaneous brain activities. As a result, the variance in averaged RSN spatial maps might be lower than the variance in time series of individual scans, which can result in higher apparent spatial correlations than temporal correlations. To control for this factor, we also calculated spatial correlations between LFP- and BOLD-derived maps for individual scans (Figs. 3B&D), and then averaged the corresponding correlation values across scans. Like scan-wise temporal correlations, scan-wise spatial correlations are significant for all LFP bands (one sample t-tests; in the M1, delta: CC = -0.37, p = 4.9 × 10− 52; theta: CC = -0.39, p = 3.9 × 10− 62; alpha: CC = -0.24, p = 5.8 × 10− 37; beta: CC = -0.15, p = 3.9 × 10− 15; gamma: CC = 0.58, p = 1.0 × 10− 56; in the ACC, delta: CC = -0.26, p = 2.6 × 10− 27; theta: CC = -0.09, p = 1.7 × 10− 6; alpha: CC = -0.04, p = 0.04; beta: CC = 0.03, p = 0.06; gamma: CC = 0.33; p = 5.4 × 10− 40). Comparisons of scan-wise spatial and temporal correlations are shown in Figs. 3E&F (paired t tests across individual scans; in the M1, delta: p = 1.01 × 10− 35; theta: p = 3.74 × 10− 50; alpha: p = 3.54 × 10− 25; beta: p = 1.18 × 10− 12; gamma: p = 7.02 × 10− 42; in the ACC, delta: p = 6.73 × 10− 21; theta: p = 5.75 × 10− 5; alpha: p = 0.74; beta: p = 7.34× 10− 5; gamma: p = 3.43 × 10− 29). These data show that even after controlling for the variance level, spatial correlations are still appreciably higher than temporal correlations.
It is likely that true LFP-BOLD temporal correlations are actually high, and lower apparent temporal correlations we observed are simply because the rsfMRI signal is noisy. Even though we cannot average rsfMRI/electrophysiology time courses across scans to reduce the noise level, we ask whether we can still estimate the true LFP-BOLD temporal correlations by assessing the noise level in our rsfMRI data. To this end, we evaluated the contrast-to-noise ratio (CNR) of our BOLD signal by leveraging the information that the spatial correlation of averaged M1 RSN maps (i.e. low-noise data) is 0.95 (Fig. 1K) while that of unaveraged RSN maps (i.e. high-noise data) is 0.58 (Fig. 3E). Based on this difference, we simulated two fixed signals with the true correlation of 0.95 (Fig. S6A). For each signal, random noise at a given CNR level was added and this process was repeated 159 times (i.e. equal to # of scans in our study). At each CNR level, the CC was calculated either based on the averaged signals from all 159 trials (i.e. low-noise data, Fig. S6A), or on signals of individual trials (i.e. high-noise data) before the resulting CCs were averaged across trials. To better mimic real BOLD/electrophysiology time series that are different across trials, for trial-wise correlations, two new signals with the same true correlation (i.e. 0.95) were simulated each trial. As expected, lower CNR gives lower apparent trial-wise CC values. Interestingly, we found that the trial-wise apparent CC of 0.58 with the true CC of 0.95 corresponds to the CNR at 1.3 (Fig. S6B-C), consistent with the CNR of BOLD contrast reported in the literature (Atkinson et al., 2008). At this CNR level, we estimated that the true BOLD-LFP temporal correlation in the M1 should be ~ 0.59 (R2 = 0.35, Fig. S6D-F) when the apparent correlation is 0.37 as measured by our data (Figs. 3A, 3E). Our simulation also indicates the difference in the number of data points (1200 in temporal correlation calculation vs 6157 in spatial correlation calculation) has a negligible influence on CC values (Fig. S6). This simulation shows that after removing the noise effect, the gamma-derived spatial map can explain up to 90% of variance in the BOLD-based seedmap, whereas the gamma power time information can only explain up to ~ 35% of variance in the BOLD time series.
All these data collectively indicate that the temporal information in the LFP signal has much lower predictive value to the local rsfMRI signal. These results are consistent with previous reports of relatively weak temporal correlations between gamma power and hemodynamic signals at rest obtained using different imaging modalities (Scholvinck et al., 2010; Winder et al., 2017).
Regressing out LFP powers does not affect RSN spatial patterns
Given the lower predictive value of the LFP power on the local rsfMRI signal, we ask to what extent the temporal information of LFP powers affect the RSN spatial patterns. The M1 (or ACC) gamma-band power, after convolving with HRF, was regressed out from the rsfMRI signals of all brain voxels. As expected, the spatial patterns of gamma power-derived correlation maps observed in Figs. 1E & 2A disappeared after the regression (Figs. 4A, 5A). However, this regression process minimally altered the M1/ACC seedmap (Figs. 4B, 5B). This result remains the same when the powers of all LFP bands were voxel-wise regressed out from rsfMRI signals using soft regression (Figs. 4C, 5C).
To examine whether the regression process is disproportionally dominated by time points with the largest LFP amplitude (i.e. outliers), we recalculated the M1/ACC seedmaps after removing rsfMRI volumes corresponding to peaks in the M1/ACC gamma power (i.e. time points with the signal amplitude above 85 percentile in HRF-convolved gamma power of each scan, Figs. 4D, 4E & 5D). The spatial similarities between the seedmaps before and after M1/ACC gamma power regression, all band power regression, or peak removal are summarized in Figs. 4F & 5E, showing that the removal of gamma power has minimal impact on the M1/ACC RSN maps.
Considering that regression is a linear process, to control for the possible nonlinear relationship between band-specific LFP powers and the rsfMRI signal, we voxel-wise calculated the mutual information between band-limited LFP powers and rsfMRI signals for all brain voxels (Fig. S7). The data show limited mutual information between any band-specific power and voxel-wise rsfMRI signals, indicating that the nonlinear component between BOLD and electrophysiological signals does not play a major role in RSN spatial patterns. These data collectively indicate the temporal fluctuations of LFP have marginal effects on BOLD-derived RSN spatial patterns.
Disparity in temporal and spatial correlations is consistent in different physiologic states
To determine whether the disparity between temporal and spatial correlations of resting-state LFP and fMRI signals we observed is a specific phenomenon under anesthesia or can be generalized to different physiologic states, we concurrently collected electrophysiology and rsfMRI data in awake rats. Despite the drastic change in the physiologic state, we again observed higher spatial correlations between LFP-derived maps (Fig. S8A) and the seedmap (Fig. S8B) in the M1, which gradually changed from negative correlations in low-frequency bands to positive correlations in high-frequency bands, revealed both in conventionally defined bands (Figs. S8C-G) and 1-Hz bands (Fig. S8H). This trend is generally maintained in the ACC, despite that the spatial correlations in low-frequency bands are diminished (Fig. S9). Similar to the results in Fig. 3, weaker but significant temporal correlations between the rsfMRI signal and HRF-convolved gamma band power were observed in unanesthetized rats (one sample t-tests; for M1, delta: CC = -0.05, p = 2.3 × 10− 4; theta: CC = -0.06, p = 4.7 × 10− 6; alpha: CC = -0.06, p = 4.6 × 10− 5; beta: CC = -0.032, p = 3.0 × 10− 3; gamma: CC = 0.04, p = 0.02; for ACC, delta: CC = 0.02, p = 0.29; theta: CC = 0.0009, p = 0.96; alpha: CC = -0.02, p = 0.31; beta: CC = 0.006, p = 0.66; gamma: CC = 0.10, p = 9.52 × 10− 7; number of scans = 50). Scan-wise comparison between temporal and spatial correlations are shown in Fig. S10. Overall, we found that the magnitudes of both spatial and temporal correlations are generally lower in the awake state compared to the light sedation state, likely due to larger variances from motion and/or other physiological fluctuations, as well as a much lower number of scans (50 awake scans, 159 light-sedation scans for M1 and 172 light-sedation scans for ACC; number of scans = 159), but the pattern of lower temporal but higher spatial correlations between gamma power and the rsfMRI signal remain consistent.
Regressing out gamma-band power or powers of all LFP bands again did not significantly alter the seedmaps for both ACC and M1 in awake rats (Fig. S11). Taken together, these results show that all major findings in lightly anesthetized rats shown above can be replicated in awake rats, suggesting that these results are not specifically related to the effects of anesthesia.
RSNs are driven by electrophysiology-invisible brain activities
Our data reveal that the LFP signal can reliably derive RSN patterns that are highly consistent with BOLD-based seedmaps. On the other hand, we found that the temporal information of the LFP signal only minimally impacts the RSN spatial pattern. To reconcile this contradiction, here we propose a theoretical model that can possibly explain this apparent paradox, described as follows:
Brain activities include components that can be measured by electrophysiology and components that are electrophysiology invisible. Electrophysiology-invisible brain activities, such as activities of nNOS neurons and/or parvalbumin (PV) interneurons (Vo et al., 2023), are actively involved in NVC, and play a dominant role in the rsfMRI signal, whereas neural activity measured by electrophysiology has minimal effects on the rsfMRI signal. At the resting state, these two components are not necessarily synchronized, leading to low temporal correlations between electrophysiology and rsfMRI signals. Another factor that can contribute to low LFP-BOLD temporal correlations is neuromodulations from distal modulator nuclei (e.g. locus coerulues and/or basal forebrain), which have strong vasoactive effects but may not commensurately affect the electrophysiology activity. Meanwhile, signaling of electrophysiology activities as well as that of electrophysiology-invisible, and thus BOLD activities are both constrained by the same anatomical pathways, allowing the two signals to separately generate similar RSN spatial patterns which can reflect both direct and indirect anatomical connectivity. The model is summarized in Figs. 6 and S13. More details of the model will be discussed in the next section.