Sleep alters neurovascular and hydrodynamic coupling in the human brain

TL;DR

This study shows that sleep changes how brain activity, blood flow, and fluid movement interact, with nonneural processes becoming more influential during sleep.

Contribution

The study reveals bidirectional coupling between neural and nonneural signals during sleep, a novel finding in human brain function.

Findings

During sleep, interactions between brain signals become more bidirectional compared to the unidirectional patterns in wakefulness.

Nonneural processes like vasomotor-driven hydrodynamic waves significantly influence brain activity during sleep.

Electrophysiological and water concentration changes predict hemodynamic changes in the awake state but not during sleep.

Abstract

Beyond its well-known effects on neuronal activity, human sleep appears also to reorganize the infraslow (<0.1 Hz) oscillation hierarchy in brain. While wakefulness is dominated by unidirectional neurovascular coupling, in which neural activity predicts hemodynamic changes, our study shows that the change in brain state from wakefulness to sleep is accompanied by increased bidirectionality in prediction patterns. Sleep-related increases in signal power were associated with bidirectional prediction patterns involving electrical activity, vascular signals, and water volume changes. These multimodal findings highlight a fundamental change in coupling dependent on brain state, suggesting that both neural and nonneural alteration contribute to sleep-related brain function.

Sleep is essential for maintaining brain tissue homeostasis, which is facilitated by enhanced cerebrospinal fluid (CSF) solute transport. Infraslow (<0.1 Hz) vasomotion, CSF flow, and electrophysiological potential all increase during sleep, but their contributions as potential drivers of CSF flow in human brain remain unknown. To investigate this, we measured these signals in healthy volunteers across sleep–wake states using functional MRI blood oxygen level-dependent (BOLD), electroencephalography, and functional near-infrared spectroscopy. We then studied the directed coupling patterns between the three signals, using phase transfer entropy. In the awake state, electrophysiological potential and water concentration changes both predicted hemodynamic BOLD changes across whole brain, reflecting classical functional hyperemia. During sleep, these interactions changed such that the net directionality was lost and the interactions became more bidirectional. Our results show that in addition to neural changes during sleep, nonneural processes such as vasomotor-driven hydrodynamic waves start to gain more impact on human brain activity.

Figures

Multimodal imaging protocol for simultaneous measurement of brain water content, electrophysiological, and hemodynamic processes. (A) Measurements were conducted in a MRI scanner using a whole-brain magnetic resonance encephalography (MREG) sequence with simultaneous direct-current electroencephalography (DC-EEG) and functional near-infrared spectroscopy (NIRS), which was specifically targeting water concentration changes. (B) Participants underwent two scanning sessions comprising functional MREG and structural scans: one during wakefulness and another during sleep. The total scan length consisted of structural scans and 10-min length fMRI scans, which were repeated during sleep scans. The pie chart illustrates the gender composition of the (n = 24) participants, while the violin plot depicts their age distribution. (C) Representative signals of infraslow (<0.1 Hz) MREG, EEG, and water NIRS from one participant along with corresponding EEG-derived sleep scoring. For visualization, NIRS signals were calibrated with respect to MREG. The pie chart shows the proportions of gathered epochs of awake, NREM-1 and NREM-2 state across all participants.

Sleep increases infraslow (<0.1 Hz) hemodynamic and electrical oscillations. (A) Time–frequency estimates illustrate the average spectral powers for magnetic resonance encephalography (MREG), electroencephalography (EEG), and functional near-infrared spectroscopy (NIRS) within sleep state-specific epochs. (B) Time-collapsed spectral estimates depict power as a function of frequency, where the infraslow frequency range (0.01 to 0.08 Hz) serves for statistical comparison. Data points corresponding to wakefulness are denoted by orange, and nonrapid eye movement sleep by light blue (NREM-1) and dark blue (NREM-2). The solid lines represent the median, while the shaded area indicates the middle 50% of data points. Statistical significance is denoted by asterisks (padj < 0.05*, 0.01, 0.001*). The EEG lacks the cardiac signal due to fMRI-related cardioballistic artifact removal.

Directed coupling patterns between infraslow brain hemodynamics, electrical, and water fluctuations. (A) Transfer entropy in the phase domain (TE) was used to calculate information transfer, utilizing the phase time-series (green) of infraslow signals (black) within sleep state-specific epochs. The sign of ΔTE indicates the main direction of information transfer, as also highlighted in the color bars. (B) Blood oxygenation level-dependent (BOLDMREG) and electroencephalography (EEG) coupling patterns were assessed in wakefulness and nonrapid eye movement (NREM)-1/2 sleep. Positive values correspond to BOLDMREG prediction of EEG, and negative values vice versa. We calculated the average TE over the whole-brain volume, where asterisks denote significance of differences (padj < 0.05*, 0.01, 0.001*). Below, the spatial distributions of statistically significant areas (padj < 0.05) are shown in comparison to awake state. (C) Similarly, we present TE between BOLDMREG and cortical water concentration changes (H2ONIRS). Here, positive values indicate BOLDMREG prediction of H2ONIRS, and negative values vice versa. We also present the whole-brain averages for BOLDMREG↔H2ONIRS coupling and spatial distribution of voxel wise differences. (D) Regional TE analysis within regions of interest (ROI) for BOLDMREG, EEG, and H2ONIRS. The box plot depicts the average TE within the ROI. (E) Topographies show average TE patterns between EEG and H2ONIRS taken over the subjects, whereas the box plot depicts average TE over the electrode space.

Spatially overlapping changes in coupling, power, and speed of blood oxygen level-dependent (BOLD) oscillations, and the hypothesized cellular mechanisms. (A) Multimodal evidence of spatial overlap in primary sensory and motor areas, showing significant increases in sleep induced infraslow (<0.1 Hz) BOLD pulsation power, wave propagation speeds (31), and changes in coupling of BOLD with both water dynamics and electrophysiological brain activity. (B) Concurrent infraslow signals, capturing BOLD signal, electroencephalography (DC-EEG), and functional near-infrared spectroscopy (fNIRS) during wakefulness and sleep. BOLD-calibrated, synchronous fNIRS signals indicate that the water (H2O) signal is anticorrelated with arterial HbO2 and BOLD signal (5, 6, 29). Changes in venous Hb concentration are less pronounced, consistent with animal studies suggesting that BOLD signal generation is driven more by speed than by volume changes in the venous compartment (32). (C) A theoretical illustration of regional electro-, hydro-, and hemodynamic interactions in wakefulness and NREM sleep. During NREM sleep, reduced orexin levels induce slow oscillations in norepinephrine (NE) levels (5, 33) and extracellular potassium [K+] concentration (34). These oscillations drive vasomotor waves detected by both BOLD and oxygenated hemoglobin (HbO2) signals (13, 14, 32), which are accompanied by anti-correlated H2O waves, as in recent animal findings (5, 10). Low NE levels reduce astrocytic volume by decreasing Na+/K+-ATPase (NAK) activity (3, 5). Pulsatile vasomotor waves cause oscillations in astrocytic endfeet and perivascular volumes (10). These oscillations facilitate hydrodynamic interstitial/cerebrospinal fluid (I/CSF) exchange via enlarged interastrocytic gaps formed by astrocytic shrinkage. This mechanism provides in theory a generation mechanism for the slow but high-voltage infraslow oscillations in EEG. Finally, the NE oscillations drive an important, increasing effect on neuronal NAK-channel activity that is opposed to the astrocytes effect closer to perivascular space. This opposing effect can explain the previously detected electro-osmotic potential difference increasing the I/CSF exchange over glia limitans explaining partially also the high-voltage DC-EEG potential (26).

Full text

Neuronal activation induces coupled arterial dilatations leading to classical activation hyperemia during the awake state (1). Sleep-associated infraslow oscillations (<0.1 Hz) in the cerebrospinal fluid (CSF) support the clearance of accumulated solutes and metabolites from the cerebrocortical parenchyma (2, 3). Norepinephrine (NE) release from cortical projections of the brainstem locus coeruleus decline in sleep and show concentration oscillations, which induce vasomotor waves in arterial wall myocytes (4, 5). The resultant slow vasomotor waves generate coupled oscillations in cerebral blood volume with inverse CSF volume changes in perivascular spaces (5?–7). The NE-driven cortical vasomotion inflates the perivascular CSF space volume and increases the thickness of astrocyte endfeet plastering the outer edge of perivascular spaces of the blood–brain barrier (BBB) (8?–10). Numerous human functional MRI (fMRI) studies have demonstrated that vasomotor wave amplitude increases with sleep depth, suggesting a direct link between sleep and vasomotor activity (7, 11????–16).

Beyond its impact on vasomotor waves, sleep also enhances infraslow oscillations in the direct current-coupled electroencephalogram (DC-EEG) (17, 18), which further modulate cortical neuronal rhythms over a wide range of frequencies (18?–20). Among these rhythms, sigma power has a notable association with the declining phase of locus coeruleus activity during transition to sleep (2, 21) and with memory consolidation after sleep (22). The infraslow EEG oscillations are further connected to changes in cerebral blood volume (CBV), extracellular pH, [K+], and permeability of the BBB (23?–25). Recent work has shown that such synchronized neuronal activity facilitates interstitial influx and efflux of electrolytes (26). Synchronous EEG & fMRI studies have also linked infraslow, synchronized neurovascular activity to enhanced solute exchange between interstitial fluid and CSF during sleep (7, 27).

A recent animal study indicated that NE oscillations coordinate vasomotor waves, CSF flow, and EEG sigma power, and that sleep alters their correlation and lag structure (5). Previous studies have often examined these factors in isolation, overlooking their potential mutual interactions. To identify the interactions underlying the sleep-associated increase in CSF flow in humans, we developed noninvasive multimodal imaging protocol not requiring tracers or contrast agents. In particular, we acquired fMRI blood oxygenation level dependent (BOLD) signals sampled at 10 Hz to avoid cardiorespiratory aliasing effects (28). Concomitant electrophysiological changes were measured using high-density DC-EEG, in conjunction with an fNIRS technology with specificity for water to assess CSF volume changes (6, 29). With this multimodal setup, we simultaneously estimated the contributions of vasomotor, electrophysiological, and CSF volume changes in the human brain.

To understand the underlying causal patterns, we studied the interactions between these three signals using information theory-based phase transfer entropy (TE) (30) across awake and sleep states. Prior studies in mice studies have demonstrated hydrodynamic modulation of blood flow oscillations (2, 5, 10). Herein, we investigated whether similar relationships are present in the awake human brain and whether these interactions are altered during sleep.

Results

To investigate the interconnected hemodynamic, water, and electrical oscillations in the human brain, we employed a multimodal neuroimaging setup (Fig. 1A). Our study included 24 healthy participants with an equal gender distribution (54% females) and a mean age of 25 y in both gender groups (Fig. 1B). Each participant completed two scanning sessions: one during wakefulness and another specifically targeting sleep. Experienced neurophysiologists (J.P., M.K.) manually classified sleep stages using the EEG recordings. We used these sleep classifications to extract in total 46 min of wakefulness, 40 min of NREM-1 sleep, and 28 min of NREM-2 sleep across all participants during active MRI scanning (Fig. 1C).

Multimodal imaging protocol for simultaneous measurement of brain water content, electrophysiological, and hemodynamic processes. (A) Measurements were conducted in a MRI scanner using a whole-brain magnetic resonance encephalography (MREG) sequence with simultaneous direct-current electroencephalography (DC-EEG) and functional near-infrared spectroscopy (NIRS), which was specifically targeting water concentration changes. (B) Participants underwent two scanning sessions comprising functional MREG and structural scans: one during wakefulness and another during sleep. The total scan length consisted of structural scans and 10-min length fMRI scans, which were repeated during sleep scans. The pie chart illustrates the gender composition of the (n = 24) participants, while the violin plot depicts their age distribution. (C) Representative signals of infraslow (<0.1 Hz) MREG, EEG, and water NIRS from one participant along with corresponding EEG-derived sleep scoring. For visualization, NIRS signals were calibrated with respect to MREG. The pie chart shows the proportions of gathered epochs of awake, NREM-1 and NREM-2 state across all participants.

Increased Infraslow (<0.1 Hz) Hemodynamic and Electrical Oscillations Occur in Sleep.

We first analyzed the signal behavior in the frequency domain to examine power levels across all modalities during sleep. Using two-min epochs of verified wakefulness and sleep (NREM-1/2), we computed the time-frequency estimates (Fig. 2A) within the 0.01 to 5 Hz frequency range. Visual inspection confirmed the presence of three distinct frequency bands: infraslow, respiratory, and cardiac ranges.

Sleep increases infraslow (<0.1 Hz) hemodynamic and electrical oscillations. (A) Time–frequency estimates illustrate the average spectral powers for magnetic resonance encephalography (MREG), electroencephalography (EEG), and functional near-infrared spectroscopy (NIRS) within sleep state-specific epochs. (B) Time-collapsed spectral estimates depict power as a function of frequency, where the infraslow frequency range (0.01 to 0.08 Hz) serves for statistical comparison. Data points corresponding to wakefulness are denoted by orange, and nonrapid eye movement sleep by light blue (NREM-1) and dark blue (NREM-2). The solid lines represent the median, while the shaded area indicates the middle 50% of data points. Statistical significance is denoted by asterisks (padj < 0.05, 0.01**, 0.001***). The EEG lacks the cardiac signal due to fMRI-related cardioballistic artifact removal.*

In the transition to sleep, there was a clear increase in the infraslow power of MREG signal (Fig. 2B), peaking at a frequency of 0.03 Hz. The group mean infraslow power rose from P_A_ = 23.3 (IQR: 23.6) during wakefulness to P_N1_ = 50.1 (IQR: 36.3) and P_N2_ = 52.2 (IQR: 37.8) during sleep (A-N1: Z_(23,20)_ = −2.59, p_adj_ < 0.05*; A-N2: Z_(23,14)_ = −3.46, p_adj_ < 0.01**). We observed a similar trend in EEG signals, with the infraslow power increasing from P_A_ = 1.49 (IQR: 2.34) during wakefulness to P_N1_ = 3.73 (IQR: 2.93) and P_N2_ = 4.48 (IQR: 3.40) during sleep (A-N1: Z_(23,20)_ = −2.42, p_adj_ < 0.05*; A-N2: Z_(23,14)_ = −2.61, p_adj_ < 0.05*). However, there were no corresponding changes in the NIRS power measurements across arousal states, with values of P_A_ = 1.15 (IQR: 1.84)*10^−11^, P_N1_ = 1.29 (IQR: 2.38)*10^−11^, and P_N2_ = 1.24 (IQR: 4.40)*10^−11^. Consistent with our expectations and earlier results, power in the infraslow range increased significantly in MREG and EEG measurements in proportion to sleep stage (13, 18).

The Main Coupling Directions between Hemodynamic, Electrical, and Water Oscillations Change during Sleep.

Based on recent observations in mice (5), we hypothesized that infraslow BOLD_MREG_ signal would be linked to electrophysiological DC-EEG and macroscopic water fluctuations H2O_NIRS_. To study the coupling patterns, we calculated the information transfer among these signals using phase transfer entropy (TE) (Fig. 3A), where TE = 0 bits indicates no connection between the phase time series.

Directed coupling patterns between infraslow brain hemodynamics, electrical, and water fluctuations. (A) Transfer entropy in the phase domain (TE) was used to calculate information transfer, utilizing the phase time-series (green) of infraslow signals (black) within sleep state-specific epochs. The sign of ΔTE indicates the main direction of information transfer, as also highlighted in the color bars. (B) Blood oxygenation level-dependent (BOLDMREG) and electroencephalography (EEG) coupling patterns were assessed in wakefulness and nonrapid eye movement (NREM)-1/2 sleep. Positive values correspond to BOLDMREG prediction of EEG, and negative values vice versa. We calculated the average TE over the whole-brain volume, where asterisks denote significance of differences (padj < 0.05, 0.01**, 0.001***). Below, the spatial distributions of statistically significant areas (padj < 0.05) are shown in comparison to awake state. (C) Similarly, we present TE between BOLDMREG and cortical water concentration changes (H2ONIRS). Here, positive values indicate BOLDMREG prediction of H2ONIRS, and negative values vice versa. We also present the whole-brain averages for BOLDMREG↔H2ONIRS coupling and spatial distribution of voxel wise differences. (D) Regional TE analysis within regions of interest (ROI) for BOLDMREG, EEG, and H2ONIRS. The box plot depicts the average TE within the ROI. (E) Topographies show average TE patterns between EEG and H2ONIRS taken over the subjects, whereas the box plot depicts average TE over the electrode space.*

We found that electrophysiological EEG changes measured during wakefulness predicted the infraslow fluctuations in BOLD_MREG_ signal (Fig. 3B and SI Appendix, Fig. S1) as expected from basic principles of neurovascular coupling. Furthermore, the water fluctuations in H2O_NIRS_ signals also predicted the BOLD signal (Fig. 3C and SI Appendix, Fig. S1, shown in blue colors) across whole brain, as indicated by a coherent negative net information transfer. In NREM1-2 sleep states, the BOLD_MREG_↔EEG and BOLD_MREG_↔H_2_O_NIRS_ coupling patterns both were altered significantly with respect to awake state (Fig. 3 B and C). The voxel level analysis showed opposing average coupling directions on cortical brain with respect to deeper subcortical brain regions.

Spatially, the area with significantly (p_adj_ < 0.05) altered BOLD_MREG_↔EEG prediction spanned a wider cortical area than the smaller but overlapping BOLD_MREG_↔H2O_NIRS_ prediction (Fig. 3 B and C). Within these significantly altered prediction regions, 81% of voxels showing increased BOLD_MREG_↔H2O_NIRS_ prediction in NREM-1 also showed a significant increase of BOLD_MREG_↔EEG, whereas there was 98% overlap for NREM-2 sleep. We further investigated if interregional delays in BOLD_MREG_ signal could explain the observed differences in cortical and subcortical prediction values. We tested this by studying voxel-wise delays within BOLD_MREG_ signal with respect to the frontal seed (SI Appendix, Fig. S2A). The mean delay maps did not correlate significantly with TE values (SI Appendix, Fig. S2 B and C). These results suggest that the observed coupling patterns did not arise from differences in these static correlation structures.

The whole brain ΔTE averages during the awake state supported the voxel level results, being clearly net negative. During sleep these averages shifted to around zero median, indicating no dominant prediction in either interaction direction. The whole brain average TE values for BOLD_MREG_↔EEG (Fig. 3B) shifted from ΔTEA = −0.11 bit (IQR: 0.07 bit) during wakefulness to ΔTEN1 = 1.810^−3^ bit (IQR = 0.09 bit) during NREM-1 and ΔTEN2 = −0.01 bit (IQR = 0.05 bit) during NREM-2. These shifts reflected statistically significant increases (A-N1: Z23,20 = −3.4, p_adj_ < 0.01**; A-N2: Z23,13 = −3.7, p_adj_ < 0.001**; N1-N2: Z20,13 = −0.06, p_adj_ = 0.96). Similarly, BOLD_MREG_↔H_2_O_NIRS_ coupling demonstrated an increasing trend over the whole-brain averages (Fig. 3C), without reaching significance. The median information transfer was ΔTEA = −0.12 bit (IQR = 0.20 bit) during wakefulness, which increased to ΔTEN1 = 0.02 bit (IQR = 0.22 bit) during NREM-1 and ΔTEN2 = 0.02 bit (IQR = 0.27 bit) during NREM-2 sleep. These cortically dominant changes were not significant for the whole brain averages: (A-N1: Z22,20 = −2.4, p_adj_ = 0.06; A-N2: Z22,13 = −1.6, p_adj_ = 0.18; N1-N2: Z20,13 = 0.4, p_adj_ = 0.87).

We hypothesized that cortical and subcortical brain regions having opposing net interaction directions in voxel level averages would cancel each other out, thus resulting in zero whole brain average TE. We further investigated the regional changes, utilizing region of interest (ROI) analysis. With BOLD_MREG_↔EEG (Fig. 3 D, Left), the net information transfer within all nine ROIs increased significantly between awake and sleep states. This effect was evident between awake-NREM1 and awake-NREM-2 states (SI Appendix, Table S1). Similarly, with BOLD_MREG_↔H2O_NIRS_ (Fig. 3 D, Right) awake state net prediction was negative within all ROIs, pointing to H2O_NIRS_ prediction of BOLD_MREG_ signal, whereas during sleep, the median centered around zero. In BOLD_MREG_↔H2O_NIRS_, there were no significant changes within the ROIs with respect to arousal state. Contrary to our expectations, the ROI analysis did not show reversed interactions directions in any of the ROIs studied during sleep, but rather indicated uniform coupling dynamics across the brain.

To further investigate the cause of the sleep-related shift in net coupling directions, we extracted the information transfer values in both interaction directions i.e., TE(x→y) and TE(y→x). By studying the separated TE values, we could then ascertain if the loss of net drive can be explained by increased information transfer in one direction, or conversely by decreased information transfer in the other direction. Our results revealed that overall information transfer in all three coupling patterns remained quite stable with respect to arousal state (SI Appendix, Fig. S3). The loss of BOLD_MREG_↔EEG net prediction direction during sleep seemed to arise from the combined effect of reduced information transfer of EEG→BOLD_MREG_ and concurrent increase in BOLD_MREG_→EEG information transfer. Similarly, loss of net prediction direction in BOLD_MREG_↔H2O_NIRS_ coupling resulted from both factors, i.e., reduced H2O_NIRS_→BOLD_MREG_ and increased BOLD_MREG_→H2O_NIRS_ information transfer.

We then focused on the final interaction pair: EEG↔H2O_NIRS_. Here, the topographic analysis demonstrated more stable coupling patterns with respect to arousal state as compared to the other interactions pairs (Fig. 3E). Generally, the central upper vertex electrodes showed continuous prediction of EEG over H2O_NIRS_, while more peripheral, distal electrodes presented an opposite prediction: H2O_NIRS_→EEG. The average ΔTE over all electrodes showed that the prediction EEG→H2O_NIRS_ was predominant, without significant alterations with respect to arousal state (Fig. 3E): (A-N1: Z22,20 = −0.14, p_adj_ = 0.9; A-N2: Z22,13 = −0.29, p_adj_ = 0.9; N1-N2: Z20,13 = 0.13, p_adj_ = 0.9).

Altogether, we found that infraslow brain hemodynamics, water volume changes, and electrophysiological changes are directionally coupled. Our results showed that the awake state is the most coherent, where water dynamics and electrophysiological brain changes both predict hemodynamic signal changes throughout entire brain. In NREM-1/2 sleep, the net prediction direction was lost due to simultaneously increased hemodynamic drive and decreased water and electrophysiological contribution.

The Sleep-Associated Increased Speed and Power of the Infraslow BOLD Overlaps with Altered Coupling Patterns in Primary Sensory Cortices.

We have previously shown with BOLD_MREG_ scans that pulsation power and speed both increase during NREM-sleep (13, 31). Here, we further explored whether these changes BOLD_MREG_ characteristics are associated with the same brain regions where the sleep-related changes in BOLD_MREG_ prediction patterns had occurred. Indeed, the increased power and speed of BOLD_MREG_ changes did overlap with voxelwise maps of the same brain regions where we observed altered coupling patterns, especially in the primary sensory cortical brain regions, posterior insula, thalamus, and upper cerebellum (Fig. 4A and SI Appendix, Fig. S4). The spatial overlap suggests either that vascular flow is less restricted in these areas or that the increased pulsation power imposes a faster flow through the porous brain tissue.

Spatially overlapping changes in coupling, power, and speed of blood oxygen level-dependent (BOLD) oscillations, and the hypothesized cellular mechanisms. (A) Multimodal evidence of spatial overlap in primary sensory and motor areas, showing significant increases in sleep induced infraslow (<0.1 Hz) BOLD pulsation power, wave propagation speeds (31), and changes in coupling of BOLD with both water dynamics and electrophysiological brain activity. (B) Concurrent infraslow signals, capturing BOLD signal, electroencephalography (DC-EEG), and functional near-infrared spectroscopy (fNIRS) during wakefulness and sleep. BOLD-calibrated, synchronous fNIRS signals indicate that the water (H2O) signal is anticorrelated with arterial HbO2 and BOLD signal (5, 6, 29). Changes in venous Hb concentration are less pronounced, consistent with animal studies suggesting that BOLD signal generation is driven more by speed than by volume changes in the venous compartment (32). (C) A theoretical illustration of regional electro-, hydro-, and hemodynamic interactions in wakefulness and NREM sleep. During NREM sleep, reduced orexin levels induce slow oscillations in norepinephrine (NE) levels (5, 33) and extracellular potassium [K+] concentration (34). These oscillations drive vasomotor waves detected by both BOLD and oxygenated hemoglobin (HbO2) signals (13, 14, 32), which are accompanied by anti-correlated H2O waves, as in recent animal findings (5, 10). Low NE levels reduce astrocytic volume by decreasing Na+/K+-ATPase (NAK) activity (3, 5). Pulsatile vasomotor waves cause oscillations in astrocytic endfeet and perivascular volumes (10). These oscillations facilitate hydrodynamic interstitial/cerebrospinal fluid (I/CSF) exchange via enlarged interastrocytic gaps formed by astrocytic shrinkage. This mechanism provides in theory a generation mechanism for the slow but high-voltage infraslow oscillations in EEG. Finally, the NE oscillations drive an important, increasing effect on neuronal NAK-channel activity that is opposed to the astrocytes effect closer to perivascular space. This opposing effect can explain the previously detected electro-osmotic potential difference increasing the I/CSF exchange over glia limitans explaining partially also the high-voltage DC-EEG potential (26).

Discussion

In this study, we used a noninvasive multimodal neuroimaging setup to investigate the brainwide coupling of infraslow (<0.1 Hz) brain oscillations of the human brain including vasomotion (BOLD_MREG_), electrophysiological signals (DC-EEG), and CSF fluctuations (fNIRS) of healthy volunteers. During the awake state, low-power infraslow electrophysiological potential and CSF oscillations both predicted vasomotor BOLD waves, consistent with neurovascular coupling (32, 35). In contrast, during sleep, the infraslow power levels increased, and the main coupling directions seen during awake state were lost. Instead, the interactions became more bidirectional, canceling out the net drive. Considering recent literature reports (2, 5), we conclude that our findings link to increased modulatory effect of vasomotor waves over CSF flow and electrophysiological oscillations in human NREM-sleep.

Hemodynamic Changes Are Coordinated by Electrophysiological and CSF Changes in Awake Human Brain.

In the awake condition, neural activation is coupled to slower increases in regional blood flow, which are traditionally thought to supply more oxygen and glucose in response to increased metabolic demand (1). However, recent findings suggest that synchronized neuronal activation following blood flow hyperemia also facilitates glymphatic solute transport and clearance, where the vasodilatory pulsations could potentially serve as an additive mechanism to accommodate increased metabolic demands (7, 26, 36?–38). Consistent with this proposition, coupling patterns in awake brain suggested sequential drive pattern in which electrophysiological activity and cortical water movement predicted BOLD changes: DC-EEG & H2O → BLOOD (Fig. 3), as has been also demonstrated in mice (5).

The water movement follows upon the neuronal [K+] release that is buffered by astrocytic K^+^IR4/5 channels, with further facilitation of passive fluid flow by neighboring astrocytic AQP4 water channels (39). The local [K+] increases are also sensed by capillary endothelia K^+^IR4/5-channels, which may facilitate upstream vasomotor dilatations via inward hyperpolarizing waves mediated by gap junctions along blood vessels (40). The occurrence of inward waves triggered by local [K+] increases offers a direct mechanism whereby NE effects on ion channels drive vascular vasomotor waves.

Reduced Coordination of Brain Hemodynamics during Human NREM-Sleep.

In sleep, the vasomotor waves were no longer driven by electrophysiological changes or by cortical water volume fluctuations. The loss of dominant direction in the coupling was associated with increased hemodynamic prediction over the electrophysiological and cortical CSF volume changes and a concurrent decrease of CSF and electrophysiological prediction over hemodynamics (SI Appendix, Fig. S3). As the information transfer baselines remained relatively constant during sleep, this implies that the interactions became more bidirectional, rather than simply decreasing the amount of overall drive. Thus, it is plausible that the direction of underlying causality oscillates over time, where the drivers take turns in providing the predominant modulation.

One potential trigger for the directionality change with time could be oscillating NE-levels. In sleep, the low level of NE release enables an increased relative activity of the Na^+^/K^+^-ATPase in astrocytes, while reducing its activity in neurons (34, 41). These opposing responses of astrocytic and neuronal Na^+^/K^+^-ATPase activity could theoretically generate an intercellular electrophysiological potential, resulting in osmotic pressure differences that would drive CSF water flow (5, 26, 34). As NE-levels slowly oscillate at ~0.02 Hz, the resulting opposite NE effects in the astrocytes and neurons could explain the increase in bidirectional effects observed during sleep (Fig. 4).

Along with oscillating NE-levels, the thickness of astrocytic endfeet at the BBB and perivascular space volume also oscillate, which could together facilitate the hydrodynamic electrolyte oscillations over the BBB glia limitans to extend further into perivascular spaces via widened interastrocytic gaps (10, 18, 20). These NE-driven vasomotor oscillatory mechanisms could underlie the increased power of DC-EEG oscillations during reduced neuronal activity in sleep found herein and as reported previously (17, 18) and could account for the increased synchrony between DC-EEG and neuronal activity during sleep (18?–20) in addition to the increased bidirectionality.

These findings, combined with recent literature reports, suggest that global vasomotor/hydrodynamic interstitial fluid oscillations under the control of NE increase during sleep, while the coupling between regional neuronal activity and hemodynamic fluctuations changes. A key factor underlying this could be marked oscillatory changes due to the bidirectional effect of global NE on neuronal and astrocyte Na^+^/K^+^-ATPase activities.

Pulsation Power and Speed Increases Overlap with Decreased Net Drive during Sleep.

Consistent with previous findings (11?–13, 15, 42), our human MREG results also showed an increase in infraslow BOLD oscillation power. Previous work has likewise shown that vasomotor waves drive the movement of injected CSF tracers within arterial wall structures (43). Furthermore, the intracranial propagation speed (31) of the vasomotor BOLD waves increases in sleep, especially in the primary sensory and motor regions, indicating an acceleration of brain water movement (Fig. 4 and SI Appendix, Fig. S4). Notably, neuronal slow delta activity also increases in the same regions, which previous animal experimentation has linked to enhanced CSF solute transport via widened interstitial spaces (13, 44). Furthermore, the electrophysiological drive of brain water tends to increase in the same parietal areas during sleep (Fig. 3E). Our results suggest that brain water increasingly flows during sleep under the drive of infraslow electrical oscillations (26).

Invasive animal experiments show that vasomotor waves are driven by infraslow fluctuations in extracellular NE levels released from projections of the locus coeruleus (2, 21), causing anticorrelated blood and CSF oscillations (5, 6). The descending phase of NE activity also correlates with the occurrence of sleep spindles and enhanced glymphatic clearance of tracers (2, 5, 21). Infraslow vasomotor waves detected in the BOLD_MREG_ (45??–48) signal also reflect reciprocal oscillations in tissue blood signal and CSF (5, 29, 45). Similarly, Hb/HbO2 levels. i.e., cerebral blood volume changes of fNIRS signal, anticorrelate with brain H_2_O volume changes (6). These multimodal lines of evidence concur in indicating that the opposing signals all arise from different compartments.

Large vasodilatory oscillations in intracranial blood volume, simultaneously captured here from BOLD signals and the arterial HbO2 signal in fNIRS (Fig. 4B), were associated with opposing water concentration changes, consistent with the Monro-Kellie doctrine, stating the sum of brain tissue, CSF, and intracerebral blood volumes must stay constant (5?–7, 10, 29).

Limitations and Future Directions.

Contrary to our expectations, the infraslow CSF oscillations measured by our H2O_NIRS_ setup did not increase during sleep. We suppose that this may be due to several methodological factors: Infrared light passes through several water compartments, which could contribute to the net water fNIRS signal. Anatomical variations across participants such as skull thickness, air sinus spaces, and CSF spaces under the frontal skull bone also directly affect the photon path length, thus altering signal sensitivity. The technical lack of multiple signal sources in fNIRS prohibit the use of averaging strategies, which might have greatly improved the signal-to-noise ratio. Due to optode placement of the fNIRS optodes, the water dynamics observed in this study are most accurately depicted in frontal brain regions, although we expect the CSF dynamics in the subarachnoid space to correlate on a larger scale.

Connectivity estimates can return spurious estimates in the presence of common drivers. In line with this, results of our previous investigation of the relationship between infraslow EEG and faster neural processing did not favor neural cortical changes as a common driver of infraslow EEG. The investigation of other possible common drivers remains a matter for future studies, where multivariate TE with increased spatial coverage of water dynamics could provide a useful framework. Future studies are needed to more thoroughly assess the effect of BOLD lag structure on signal prediction values. Furthermore, the use of longer data epochs may allow for the detection of more nuanced coupling effects and enable the examination of temporal changes in signal prediction patterns.

Materials and Methods

Experimental setup. This study was approved by Regional Ethics Committee of the Northern Ostrobothnia Hospital District. We obtained written informed consent from all participants in accordance with the Declaration of Helsinki. The datasets utilized in this study were originally gathered as part of our prior investigation (13, 14, 18) involving 24 healthy controls (13 females, 11 males). Participants underwent scanning twice across separate sessions: once during wakefulness and once during sleep. The total duration of awake scans was approximately 30 min, while the sleep scans lasted for about an hour. The measurement protocol included structural sequences and fMRI MREG sequences lasting 10 min each, which were repeated during sleep scans. To increase sleep pressure and enabling faster onset of the sleep recordings, 13 subjects underwent one night of sleep deprivation, which was monitored with smart rings (Oura Health Oy). For each participant, 2-min data epochs were extracted, as these represented the longest continuous segments available across subjects and sleep stages. Data outside these epochs were excluded to avoid unequal weighting across participants. Our study employed a multimodal neuroimaging setup (49) including an MRI scanner to capture functional ultrafast MREG and structural image series. Spontaneous brain activity was recorded with EEG, which was further used to derive sleep state–specific information. Sleep scoring was performed in 30-s epochs by experienced neurophysiologists (J.P., M.K.) according to AASM criteria. Functional NIRS was used to measure changes in macroscopic water concentrations. All modalities were synchronized with the MRI scanner’s optical timing pulse. Calculations were performed using MATLAB (v.R2023b, MathWorks).

MRI Acquisition and Preprocessing.

Functional and structural images were acquired using a Siemens MAGNETOM Skyra 3T scanner equipped with a 32-channel head coil. For structural 3D MPRAGE scanning (TR = 1,900 ms, TE = 2.49 ms, TI = 900 ms, FA = 9°, FOV = 240 mm), we used a slice thickness of 0.9 mm. Functional imaging was performed using an ultrafast MREG-sequence (TR = 100 ms, TE = 36 ms, FA = 5°, FOV = 192 mm), which employs k-space undersampling to reach a sampling frequency of 10 Hz with voxel size of 3 mm (50). We set a crusher gradient to 0.1, which was optimized for detecting physiological signal sources, while mitigating slow drifts and stimulated echoes.

Reconstruction of MREG images involved the utilization of L2-Tikhonov regularization, wherein a lambda value of 0.1 was determined using the L-curve method (51). Dynamic off-resonance correction in k-space was implemented to reduce B0-field artifacts and mitigate respiratory motion. Image preprocessing followed the standardized FSL (Functional MRI of the Brain’s software library) preprocessing pipeline (52). A high-pass filter was used with a cut-off frequency of 0.008 Hz. Subsequently, we performed motion correction followed by brain extraction. To eliminate artifactual spikes, datasets were despiked (53). The structural 3D MPRAGE images were used in the registration of functional datasets into the MNI152 standard space (Montreal Neurological Institute). Continuous 2-min segments of wakefulness and NREM-1/2 sleep were identified using EEG-derived sleep scores (14).

EEG Acquisition and Preprocessing.

EEG recordings were acquired using a GES 400 (Magstim EGI) system, which included a direct current-coupled amplifier (Net Amps 400), and a high-density 256-electrode system (HydroCel Geodesic Sensor MR net) with the electrode “Cz” serving as the reference channel. The recordings were conducted at a sampling rate of 1 kHz, except for three sleep and five awake subjects in whom the experimentalists had inadvertently selected 250 Hz. Prior to measurements, we performed a visual inspection of signal quality and electrode impedances.

The initial steps in data processing involved the removal of gradient and ballistocardiographic artifacts through template subtraction (54, 55); Brain Vision Analyzer v.2.1, Brain Products). Subsequently, bad channels were excluded based on the following criteria: SD exceeding 2,000 µV, average correlation with neighboring electrodes falling below 0.1, or electrode impedance surpassing 1 MΩ. Upon identifying bad channels, spherical interpolation was used to replace the removed channels (56). EEG signals were decimated to 10 Hz to align with the sampling frequency of other modalities.

NIRS Acquisition and Preprocessing.

Our functional NIRS device (29) utilized a frequency-coding technique in which the emitted light was modulated at specific frequencies for each wavelength. High-power LEDs generated monochromatic light at wavelengths of 690, 830, and 980 nm, facilitating measurements of Hb, HbO, and H_2_O, respectively. Specifically, we selected the 980 nm wavelength for its high absorbance of water, ensuring high sensitivity in water dynamics measurement. At the receiver optode, the light was demodulated again with the corresponding frequencies. To maintain consistent optode positioning, we separated the source and receiver optodes on both sides of the EEG electrode “Fp1” at a fixed distance of 3 cm, thereby ensuring a constant distance between measurements, while enabling near-infrared light to reach the cerebral cortex (29).

Initially, NIRS transmittance values were used to compute HbO, HbR, and water concentrations following the modified Beer–Lambert law (29, 57). NIRS concentration signals were decimated to 10 Hz to match the sampling rates of all three modalities.

Estimation of Spectral Properties.

We employed complex Morlet wavelets in wavelet convolution to conduct time–frequency spectral analysis (Fig. 2 A and B). These wavelets consist of sine waves at various frequencies tapered by a Gaussian window, thereby providing temporal specificity. Convolving time series with each frequency’s wavelet transforms the data into frequency-specific power by computing the squared magnitude of the convolution results. Unlike rectangular windowing, the Gaussian window reduces ripple effects typically associated with sharp edges in the kernel signal. We chose wavelet convolution due to its ability to focus on time-domain changes, its computational efficiency, and its adherence to the assumption of stationarity, contrary to conventional spectral methods (58). We applied a logarithmic frequency range extending from 0.01 Hz to 5 Hz with 50 steps, having selected the 5 Hz upper limit as it represents the Nyquist frequency of the MREG recordings. To mitigate edge effects that could potentially contaminate the time series, particularly in low-frequency filtering, we implemented the mirroring technique. The number of wavelet cycles (N = 8) was held constant, regulating the trade-off between frequency and temporal precision. Given that we had nearly 70,000 spatially correlating signals in each MREG dataset, we employed a simple random sampling technique without replacement to reduce computational cost. Here, we used 5% (3,405 voxels) of all voxels, ensuring that each brain voxel had an equal probability of being chosen. Finally, we calculated power in the frequency range of 0.01 to 0.08 Hz using integral approximation and rectangular windowing. We averaged the spectral power estimates over the MREG voxels and EEG electrodes to obtain global spectra.

Inferring Directed Coupling Patterns.

Directed coupling patterns were analyzed in the infraslow frequency range, using the 2-min sleep state–specific epochs (Fig. 3). With sampling rates matched across modalities, we constructed a finite impulse response (FIR) filter with Hamming windowing, filter order 3,000, and cutoff frequencies 0.01 and 0.08 Hz to bandpass filter the signal into the infraslow frequencies. The chosen frequency band allowed the inclusion of slow hemodynamic processes, but leaving out potential slow respiration effects and faster cardiac effects. The mirroring method was used to ensure that no edge effects remain after temporal filtering, and zero-phase recursive filtering was used to avoid phase offset and distortions.

Standard correlative metrics can fail to inferaccurately the directed interactions, especially in the presence of nonlinear interactions. To address this limitation, we used phase transfer entropy (TE), a nonlinear extension of the classical Granger causality, to quantify putative causal interactions among physiological signals (30). TE is a dynamic and directed measure of information transfer (30, 59). Essentially, TE represents the average information from a source signal that aids in predicting the next value of the target concerning its past. We calculated TE in the phase domain, which has been shown to be more robust in comparison to real valued TE in the presence of environmental noise (30). We applied Hilbert transform for the infraslow signals to compute the analytical signals $[eqn]$ . From these analytical signals, we extracted the instantaneous phases as $[eqn]$ .

In calculating TE, we chose a discrete estimator to represent the state space, being the simplest form of estimators. Marginal entropies, joint entropies, and phase transfer entropies were defined as

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where t’ refers to past time points t’ = t-δ, and δ is the analysis lag.

If the source–target delay is set to maximize information transfer, it aligns with the actual causal delay under simple conditions. However, predictive information can be accurately estimated over a wide range of lags (30). Here, we assumed a constant delay of one cycle, in consideration of the computational cost of identifying the maximal information transfer. We calculated pairwise TE for each MREG voxel, including all EEG electrodes and water concentration signals. TE values were then averaged over the EEG electrode space and for whole brain averages also over the voxel space, aiming to reduce data dimensionality. We used FSLs MNI structural atlas to divide the brain into 9 bilateral ROIs: caudate, cerebellum, frontal lobe, insula, occipital lobe, parietal lobe, putamen, temporal lobe, and thalamus. We then extracted the regional ΔTE metrics and averaged over the electrodes and voxels coinciding within the ROIs. To infer the net direction of the information transfer, we took the difference as $[eqn]$ .

Vasomotor Wave Propagation Speed.

To investigate vasomotor wave propagation within the human brain, we applied optical flow analysis to the MREG data, focusing specifically on infraslow oscillatory components associated with vasomotor activity. Optical flow analysis has been previously used for tracking motion within imaging data for different physiological pulsations including cardiac (60) and respiratory pulsations (61). By leveraging this method, we aimed to capture the spatiotemporal dynamics of vasomotor waves (31). We quantified the velocity and spatial extent of vasomotor wave propagation, thereby characterizing how these waves propagate across cortical and subcortical regions.

Pulse-triggering was guided by identifying regions demonstrating the highest correlations with vasomotor activity, following established methodologies. Specifically, the posterior cingulate cortex emerged as the region most strongly associated with vasomotor oscillations, leading to its selection as the primary reference region for triggering vasomotor wave analysis. This approach provided a precise framework for delineating the onset and propagation of vasomotor oscillations within the brain’s functional architecture.

Statistical Analysis.

We had not made a priori statistical power calculation to predetermine the sample size for this study but used the material at hand. The alpha level was set to 0.05 for all statistical tests. We investigated whether infraslow power differed between epochs of wakefulness and the combined sleep states (NREM-1/2). A Wilcoxon rank-sum test (two-tailed) was employed, with the null hypothesis being that awake and sleep epochs had equal infraslow power (Fig. 2B). The false discovery rate (FDR) was controlled using FDR correction (62) to account for multiple comparisons. To compare voxelwise group differences in TE, we conducted a randomization test with 5,000 permutations, shuffling the arousal state labels for awake; NREM-1, awake; NREM-2, and NREM-1; NREM-2 (Fig. 3 B and C). Here, we used the threshold-free cluster enhancement method to address multiple comparisons. In comparison of whole-brain average TE and within ROI TE, we used a Wilcoxon rank-sum test (two-tailed), with the null hypothesis being that arousal states had equal median information transfer (Fig. 3 B–E). Here as well, we applied correction for the FDR. For comparison of BOLD power and speed maps between awake and sleep states, we used a randomization test along with threshold-free cluster enhancement (31).

Supplementary Material

Appendix 01 (PDF)