Changes in cross-frequency coupling following closed-loop auditory stimulation in non-rapid eye movement sleep.

Regional changes of non-rapid eye movement (NREM) sleep delta and sigma activity, and their temporal coupling have been related to experience-dependent plastic changes during previous wakefulness. These sleep-specific rhythms seem to be important for brain recovery and memory consolidation. Recently, it was demonstrated that by targeting slow waves in a particular region at a specific phase with closed-loop auditory stimulation, it is possible to locally manipulate slow-wave activity and interact with training-induced neuroplastic changes. In our study, we tested whether closed-loop auditory stimulation targeting the up-phase of slow waves might not only interact with the main sleep rhythms but also with their coupling within the circumscribed region. We demonstrate that while closed-loop auditory stimulation globally enhances delta, theta and sigma power, changes in cross-frequency coupling of these oscillations were more spatially restricted. Importantly, a significant increase in delta-sigma coupling was observed over the right parietal area, located directly posterior to the target electrode. These findings suggest that closed-loop auditory stimulation locally modulates coupling between delta phase and sigma power in a targeted region, which could be used to manipulate sleep-dependent neuroplasticity within the brain network of interest.

Introduction

A substantial body of evidence highlights the importance of non-rapid eye movement (NREM) sleep for experience-dependent neuroplasticity and associated improvements in behavioral performance after a bout of sleep[1,2]. NREM sleep electroencephalographic (EEG) recordings are hallmarked by neocortical slow waves (1–4 Hz) and sleep spindles, the waxing and waning activity in the sigma frequency band (10–16 Hz). It is well established that slow waves orchestrate faster oscillations, with the up-phase of slow waves being associated with increased sigma band activity, and sigma oscillations, in turn, group gamma activity in their troughs[3,4]. This precise cross-frequency phase-amplitude coupling is thought to organize communication across multiple brain areas and to facilitate coordinated information processing in brain networks during NREM sleep[5].

Not only the magnitude and spatial distribution of slow waves and sigma activity[3,6-10], but especially their temporal coupling seems to contribute to sleep-dependent benefits of behavioral performance[11,12]. Similar to phase-amplitude coupling between slow waves and sigma, modulation of theta band activity (5–8 Hz) by the phase of slow waves has been described in several studies. Both in surface EEG and in intracranial (iEEG) recordings, theta power was increased during the transition to the down phase of slow waves, with the last trough of theta-burst coinciding with the slow-wave down-phase[5,13-15].

Moreover, theta band activity in NREM sleep may play a role in sleep-dependent memory consolidation in humans. NREM theta oscillations were associated with the replay of the wake neuronal activity patterns during slow-wave sleep[16,17]. Furthermore, theta-range oscillations were linked to the benefits of targeted memory reactivation, a technique where learning-related cue stimuli (e.g. an odor or a sound) are presented during NREM sleep to promote sleep-dependent memory consolidation[17-21]. Thus, clear evidence exists both for coupling between slow waves and NREM theta and the association of theta band activity with memory replay. Nevertheless, theta is still not included as a primary component in the theoretical framework describing the functional importance of the temporal interaction between slow waves and faster rhythms for sleep-dependent neuroplasticity. A thorough description of the coupling between slow waves, theta, and spindles within the same data set might provide novel insights into the role of theta activity in these processes.

The manipulation of slow-wave activity during sleep by various techniques was instrumental in establishing causal relationships between sleep and learning-related processes[22-26]. However, whether cross-frequency coupling is critical for such causal relationships between sleep and learning is unknown. A first step would be to show that it is possible to non-invasively manipulate cross-frequency coupling. Closed-loop auditory stimulation, a novel approach to modulate slow waves, can be used to collect such evidence. It has been shown that these brief auditory stimuli, time-locked to the up phase of endogenous slow waves, not only induced trains of high-amplitude slow waves but also boosted sigma and theta activity[27-31]. The analyses of evoked event-related potentials (ERP) and event-related spectral perturbations (ERSP) show that they share the primary characteristics of spontaneously generated slow waves. Specifically, these evoked slow waves follow a similar time–frequency pattern, where the transition to a down-phase is associated with theta activity (~ 0.4 s after stimulus onset) and the up-phase coincides with a peak of sigma activity (~ 1 s after stimulus onset). Interestingly, this theta and sigma power enhancement coupled to evoked slow waves is not only found in studies using closed-loop auditory stimulation with pink noise bursts[29,31,32], but also in studies using targeted memory reactivation paradigms[19,33]. The functional relevance of this similarity between a specific reactivation of activity (by cueing) and the seemingly unspecific boost of slow waves by closed-loop auditory stimulation is unknown.

One important aspect to consider in this context is the local regulation of slow-wave activity. For example, sleep-dependent plasticity is reflected in locally enhanced slow-wave activity during sleep following intensive training on a particular task during wakefulness[34,35]. This observation suggests that the “need for sleep” increases as a function of learning-related brain plasticity. More recently, similar local effects of pre-sleep learning were also found to be related to cross-frequency coupling between slow waves and sigma activity[36-38]. Hence, it might be that local enhancement of power and cross-frequency interaction of the main sleep EEG rhythms in the relevant brain regions is more crucial for the post-sleep performance gains as compared to global changes in these parameters. The critical question that arises is whether closed-loop auditory stimulation during sleep results in global or local effects.

Thus, the goal of our analysis was to test whether slow-waves, theta and sigma oscillations and their interactions can be modulated specifically in a particular brain region by closed-loop auditory stimulation. Spatially restricted modulation of cross-frequency coupling could be potentially used as a tool to interact with sleep-dependent memory formation within the brain network of interest.

Results

We analyzed high-density EEG recordings of 9 subjects (18 full-night recordings). Subjects participated in two-night sessions separated by 1 week: non-stimulation (SHAM) and closed-loop auditory stimulation in the up-phase (STIM) of real-time detected slow waves. Real-time slow-wave detection during sleep stages N2 and N3 was performed as described in our previous publication[23]. In the STIM condition, auditory stimuli (50 ms bursts of 1/f pink noise, inter-stimulus interval ≥ 2 s) were played whenever the EEG signal in the C4 electrode (the “target channel”) crossed a default threshold of + 30 µV (Fig. 1, see Materials and Methods). Subjects were blind to the order of experimental nights and objective sleep parameters did not differ between STIM and SHAM nights except for the total time in bed (Supplementary Table 1). Importantly, there was no difference in the number of arousals (t(8) = 0.97, p = 0.3) during the stimulation period in the STIM condition (9 ± 3) and the corresponding period of the SHAM condition (11 ± 4).

Figure 1

Closed-loop auditory stimulation protocol. During the STIM experimental nights, acoustic stimulation (50 ms bursts of 1/f pink noise, inter-stimuli interval ≥ 2 s) was delivered whenever the EEG signal crossed a threshold of  + 30 µV (dashed lines) in the electrode C4 electrode located over the right sensorimotor area (red electrode on the high-density EEG cap).

ERP and power analysis

We first explored the ERP and time–frequency representations (TFR) following the auditory stimuli presentation in the STIM condition and contrasted them to the SHAM condition, averaging across all channels. As shown in Fig. 2, the presentation of auditory stimuli was followed by a pronounced ERP response, resembling a classical K-complex. The time-locked response contained all of the classical ERP components, including the P200, N350, N550 and P900[39]. When comparing these ERPs with spontaneous slow waves in SHAM, we observed a significant change in the latency of the down-phase (t(8) = −4.75, p = 0.001), occurring earlier in the STIM condition (0.55 ± 0.04 s, M ± STD) compared to the SHAM condition (0.73 ± 0.09 s). A change in latency was also observed for the up-phase (t(8) = −6.03, p < 0.001): again, the up-phase occurred earlier in the STIM (1.1 ± 0.07 s) than in the SHAM condition (1.50 ± 0.16 s). Not only the latency of down- and up-phase was shortened, but also the transition between down- and up-phase occurred faster (t(8) = −3.44, p = 0.009), with a mean inter-peak latency of 0.55 ± 0.06 s in the STIM condition and 0.76 ± 0.16 s in the SHAM condition.

Figure 2

Total power changes following the stimulus presentation during auditory closed-loop stimulation. Superimposed grand-averaged ERP waveforms (in C4 target channel) for STIM (solid line) and SHAM (dashed line) conditions in an overlay, with (left) time–frequency power plots (relative change in power between STIM and SHAM collapsed across all channels) and (right) time–frequency t-value plots (shaded area indicates non-significant difference between conditions, cluster corrected two-sided p < 0.05). Note that the stimulus presentation was associated with a power increase in the theta band during the transition to the down phase of slow waves, followed by a power increase in fast sigma and gamma bands during the up-phase of slow waves.

In the next step, we assessed total power differences between STIM and SHAM conditions across all frequencies from 0.75 to 40 Hz in the 2 s after stimulus onset using cluster-based permutation tests (Fig. 2, TFR plots, average of all channels). This analysis showed that the activity in the delta (~ 0.5–2 Hz) and theta range (~ 3–10 Hz, prior to the slow-wave down-phase) was enhanced in STIM compared to SHAM by 48 ± 35% and formed a cluster (p = 0.02, d = 1.31). Similar to previous findings, our TFR results have shown that up-phases of slow waves were associated with a robust increase of 44 ± 26% in fast sigma power (p < 0.001, d = 0.77), with a center frequency of ~ 16 Hz and extending to higher frequencies (slow gamma, ~ 30 Hz). Interestingly, at the end of the 2-s analysis window, there was also a suppression of activity in the sigma band following the enhancement (p = 0.03, d = −0.38; mean decrease −13 ± 6%), which is in line with previous findings on refractory periods of spindle occurrence[40].

To determine the spatial distribution of the observed effect, we characterized stimulation-related topographical changes in different frequency bands. A significant power enhancement within the delta band (0.75–3 Hz) was observed in frontocentral areas during the first second after stimulus onset (p = 0.02, d = 1.05, mean increase 41 ± 34%). The widespread stimulation-related theta (4–8 Hz; p = 0.003, d = 0.98, mean increase 40 ± 19%) and sigma (13–17 Hz; p < 0.001, d = 0.81, mean increase 51 ± 34%) power enhancement also involved frontal and central areas. The time-resolved topographies in Fig. 3 show that in both frequency bands the power increase arises as a local cluster of central channels, which spreads over frontal and parietal areas and lasts the longest over frontal areas. The cluster of power enhancement was distributed differently for theta and sigma components: in the theta band, the focus was centered around the vertex while in the sigma band the focus was in posterior channels over the parietal lobes.

Figure 3

Topographical distribution of delta, theta and sigma power changes between STIM and SHAM conditions. Power value ratios (STIM/SHAM) were averaged in 200 ms time windows, starting from 200 ms after stimulus onset (start of the stimulation effect in the delta-theta range according to Fig. 2). White dots indicate significant changes, cluster corrected two-sided p < 0.05. The red dot indicates the target channel C4.

In summary, this analysis shows a generalized enhancement of power in delta, theta and sigma bands following closed-loop acoustic stimulation, distributed predominantly in frontal and central areas.

Cross-frequency coupling analysis

To characterize delta-theta and delta-sigma coupling in STIM and SHAM conditions, we first compared the preferred delta phase and modulation index (MI) in both conditions in 9 representative channels across the cortex (F3, Fz, F4, C3, Cz, C4, P3, Pz, P4).

Coupling strength

To determine if delta-theta and delta-sigma coupling strength was affected by up-phase stimulation, we computed the MI at each channel in both conditions. The MI provides information only about the degree to which the amplitude of theta or sigma is modulated by the phase of delta oscillations. We normalized the coupling strength using a permutation-based reshuffling approach, resulting in a z-scored coupling strength (zMI), independent from differences in absolute power between conditions, channels, and individuals.

In the 9 reperesentative channels, the coupling strength between delta and theta was significantly greater than zero both in STIM (all t(8) > 5.32, padj < 0.001, FDR corrected) and SHAM (all t(8) > 5.66, padj < 0.001). The coupling strength between delta and sigma was also significantly greater than zero in STIM (all t(8) > 2.97, padj < 0.02) and in SHAM (all t(8) > 2.23, padj < 0.05). The effect sizes were slightly larger in frontal channels for theta coupling (STIM, frontal channels: d = 3.17, central: d = 2.45, parietal: d = 2.12) and in central channels for sigma coupling (STIM, frontal channels: d = 1.14, central: d = 1.89, parietal: d = 1.07). Similar patterns of interregional differences were observed for the SHAM condition. These results were confirmed later by the topographical distribution of zMI (Fig. 5).

Figure 5

Topographical distribution of coupling strength. Delta-theta and delta-sigma coupling (z-scored modulation index, zMI) topography for STIM (left) and SHAM (middle) conditions and the ratio of STIM to SHAM (right). White dots indicate significant changes (cluster corrected two-sided p < 0.05), the red dot indicates the target channel C4.

Comodulograms in Fig. 4 demonstrated region- and frequency-specific effects of the stimulation on coupling strength. A significant decrease in delta-theta coupling was observed in frontal channels F3 (p = 0.004, d = −1.62) and F4 (p = 0.01, d = −1.45) in the STIM condition as compared to SHAM. In the STIM condition, coupling strength for delta and sigma bands was increased in the parietal electrode P4 (p = 0.004, d = 1.2), located posterior to the target channel C4. Interestingly, in the vertex region, there was a significant decrease in coupling strength between the upper delta (~ 3–4 Hz) and sigma (Cz, p < 0.001, d = −1.85) bands.

Figure 4

The z-scored modulation index as a function of amplitude (3–40 Hz) and phase (0.75–5 Hz) for nine electrodes. Non-shaded areas indicate statistically significant changes in cross-frequency coupling between STIM and SHAM conditions (cluster corrected two-sided p < 0.05). Reduction in coupling strength for delta and theta was observed in frontal channels F3 and F4 while coupling strength for delta and sigma bands was increased in the parietal electrode P4, located posterior to the target channel C4.

Topographical distribution of zMI changes

Topographical distributions of delta-theta and delta-sigma coupling in STIM and SHAM conditions are shown in Fig. 5. Delta-theta coupling was overall stronger in frontal sites both in STIM and SHAM conditions. For delta-sigma coupling, the focus was located over parietal sites. These results are in agreement with interregional differences in the effect sizes observed in 9 representative channels above.

We then compared zMI topographies for the two conditions. The up-phase stimulation caused a significant decrease of delta-theta coupling by 31 ± 18% in frontal channels (p = 0.015, d = −1.64), the area of the strongest default coupling both in STIM and SHAM conditions. For delta-sigma coupling, we observed a significant increase by 124 ± 105% over the right parietal area, located directly posterior to the target channel (p = 0.01, d = 1.15).

Contribution of the phase-locked response to changes in delta-theta coupling

The overlay of TFR with ERP responses in Fig. 2 indicates that the theta increase observed prior to the down-phase in STIM might contain non-phase-locked power changes (induced responses), as well as phase-locked power changes (evoked responses, including P200 and N350 components of the ERP). The majority of the theta changes that are observed 500 ms after stimulation are non-phase-locked, as can be seen on the induced power TFR in Figure S1 (for more details on phase-locked and non-phase-locked power analysis see Frequency analysis of the Methods section). A similar assumption can be made about activity in the delta range, containing both ongoing non-phase-locked activity overlaid with phase-locked stimulation-related components N550 and P900.

Along these lines, the observed decrease in coupling between delta and theta in frontal regions might be explained by an altered relationship between these two frequency bands due to the presence of evoked delta and theta responses. To test if the phase-locked activity in delta or theta bands is contributing to the observed decrease in coupling, in addition to zMI computed using total power (1), we calculated zMI using only: (2) broadband induced power (by subtracting the broadband ERP), (3) delta induced power (by subtracting the delta-filtered ERP) or (4) theta induced activity (by subtracting the theta-filtered ERP). Accordingly, zMI computed from total power includes both evoked delta and theta, whereas zMI computed from induced broadband power would not contain either of them.

A three-way ANOVA with factors Condition (STIM vs. SHAM), EvokedDelta (present or subtracted) and EvokedTheta (present or subtracted) yielded a significant interaction between Condition and EvokedDelta [F(1,8) = 17.96, p = 0.003, η2 = 0.69]. The effect of EvokedDelta was significant for STIM [F(1,8) = 32.76, p < 0.001, η2 = 0.8], but not for SHAM [F(1,8) = 1.16, p > 0.3, η2 = 0.13]. Thus, only in the STIM condition, stronger coupling was observed when evoked delta activity was subtracted from the data (73.47 ± 16.24 zMI) as compared to the coupling calculated including the evoked delta response (58.96 ± 17.54 zMI). The analysis also revealed a significant main effect of Condition on coupling strength [F(1,8) = 9.74, p = 0.014, η2 = 0.55], resulting from an overall weaker coupling in the STIM condition as compared to SHAM. The main effect of EvokedDelta was also significant [F(1,8) = 18.50, p = 0.003, η2 = 0.7]. The main effect of EvokedTheta was not significant, indicating that subtraction of the evoked theta did not affect zMI. No other significant interactions between factors were observed. Taken together, the ANOVA analysis demonstrated that the stimulation-related phase-locked delta could mediate the frontal decrease of cross-frequency coupling observed after up-phase stimulation.

Delta phase preference

Theta and sigma oscillations were preferentially expressed during different phases of delta waves. Theta activity was maximal in the transition to the down phase of delta waves both for STIM and SHAM (average phase across 9 electrodes, averaged across subjects: 119° ± 10° [circular mean ± SD] in STIM and 119° ± 7° in SHAM; Figure S2). On the other hand, for the sigma band, the activity was preferentially expressed close to the up-phase of delta waves (−19° ± 32° in STIM and −5° ± 21° in SHAM; Figure S1). The phase distribution averaged across channels was non-uniform for theta and sigma power in both experimental conditions (Rayleigh test: all p < 0.01).

The circular distributions in Figure S2 indicate the spatial extent of the preferred delta phase for theta and sigma frequencies. Because both theta and sigma maximal power values were scattered within half a cycle of the delta oscillation, we used a parametric two-way ANOVA with factors Location (frontal, central, posterior) and Condition (STIM vs. SHAM). Significant topographical differences in the preferred delta phase were observed both for theta and sigma. The distribution of the preferred phase formed a frontoparietal gradient, with the power maximum occurring earlier in the cycle in parietal areas. Thus, the maximum of theta power was located closer to the down phase in frontal sites and closer to the up-to-down transition in parietal sites [main effect of Location, F(1.56,12.21) = 46.34; p < 0.001, η2 = 0.85; preferred delta phase: 143° ± 3° in frontal, 118° ± 3° in central and 93° ± 5° in parietal sites]. Similar topographical differences were observed for sigma: the maximum of sigma power was located slightly after the up-phase of delta waves in frontal regions and slightly before in parietal sites [main effect of Location, F(1.48,11.86) = 22.02; p < 0.001, η2 = 0.73; preferred delta phase: 28° ± 13° in frontal, −13° ± 8° in central and −30° ± 5° in parietal sites]. There was no difference in the mean coupling direction between STIM and SHAM condition for the theta band [main effect of Condition, F(1,8) = 1.61; p > 0.2, η2 = 0.17]. However, across all locations, the maximum of sigma power occurred earlier in STIM (−12° ± 7°) compared to SHAM (2° ± 9°), as indicated by the significant main effect of Condition [F(1,8) = 8.29; p = 0.021, η2 = 0.51]. No significant interaction between the factors Condition and Location was observed. In sum, acoustic stimulation significantly affected the preferred delta phase for sigma, but not for theta. In addition, we observed changes in phase preference across 9 representative electrodes.

Discussion

Although the majority of closed-loop auditory stimulation studies are still conducted in lab environments, recently, the first systems for unsupervised closed-loop enhancement of slow-wave activity have been developed[41,42] with the goal to boost memory consolidation for prolonged periods and in clinical populations. Given the importance of local sleep characteristics for use-dependent neuroplasticity and considering the translational trend in the auditory closed-loop stimulation field, it is important to investigate whether slow waves, theta and sigma oscillations and their interaction can be manipulated selectively in a particular brain region. Here, we demonstrate that while closed-loop auditory stimulation globally enhances power in delta, theta and sigma bands, changes in cross-frequency coupling between these oscillations are expressed more locally. In particular, a significant increase in coupling between slow waves and sigma power was observed over the right parietal area, located directly posterior to the target electrode. This demonstrates that a local modulation of the coupling between slow waves and sigma band activity by closed-loop auditory stimulation could be used to manipulate sleep-dependent memory formation within the brain network of interest.

Similarities between time–frequency profiles of stimulation-related responses and endogenous slow waves

By applying auditory stimulation with an inter-stimulus interval of ≥ 2 s, we had a sufficiently long-time window to perform a detailed analysis of the ERP responses and associated spectral changes. Evoked ERP responses contained all the time-locked spectral components typical for stimulation-evoked K-complexes (P200, N350, N550, and P900)[43]. A time-resolved topographical analysis of the stimulation-related delta changes demonstrated a similar smooth transition between the near-vertex N350-like distribution and the frontocentral N550-like distribution within the first second after stimulation as found previously[44]. K-complexes and the recently proposed type I slow waves[45] are thought to be generated in sensorimotor areas and the medial parietal cortex and involve predominantly frontocentral regions[43, 45–47]. Interestingly, stimulation-induced delta changes observed in our study were also distributed mainly in central areas, indirectly suggesting recruitment of the same neural mechanisms of slow-wave generation. Thus, a plausible interpretation of the commonly found closed-loop stimulation effects is that by applying stimulation time-locked to a particular phase of slow waves, further synchronization of endogenous slow waves is promoted by means of the externally evoked K-complex[28]. This conclusion is in agreement with the interpretation of closed-loop auditory stimulation effects initially proposed by Bellesi et al. in 2014[48]. However, we cannot exclude that the frontocentral predominance of the observed changes in delta power may be related to the location of our target electrode C4. Shifting the target electrode to e.g. the occipital cortex would be an interesting probe of this idea for future experiments.

Similar to other studies using auditory stimulation during sleep[19,29,31-33], evoked responses within the delta range were associated with non-phase-locked activity in other frequency bands, namely theta and sigma. Changes in theta and sigma power were broadly distributed over the scalp and there was no target-electrode-specific local effect. However, there was a clear temporal relationship between slow waves and power in the theta and sigma frequency range as previously described in the literature and found not only after stimulus presentation, but also during spontaneous slow waves[13,14,32,36,37,49,50]. Interestingly, a very similar time–frequency pattern was observed in a TMS study, where single pulses over the primary motor area boosted sigma and theta power and these changes had a similar temporal relationship with down-phases (tms-N400) and up-phases (tms-P1000) of the evoked slow waves[51]. Thus, theta-sigma enhancement associated with evoked K-complex-like responses might be of particular interest for understanding mechanisms underlying closed-loop auditory stimulation effects, as it represents the general “default” mechanism of brain responses to external stimulation.

Effects of closed-loop auditory stimulation on coupling between the delta band and faster oscillations

Most recently, it has been shown that not only individual properties of slow waves and spindles can be associated with benefits of sleep for recovery and memory formation, but also the interaction between them[6,37,52,53]. The relationships between slow waves and sigma are well described in the literature. It is widely accepted that cortical down-states via cortico-thalamic connections are able to trigger activity in the sigma range which in turn is associated with the up-phase of delta waves[4]. Notably, there are very few scalp-EEG studies describing the generative mechanism of theta oscillations in NREM sleep and their relationships with slow waves.

In our study, we investigated the cross-frequency coupling between waves in the delta frequency range and faster oscillations in theta and sigma bands. Overall, despite the absence of local effects in power changes, auditory stimulation clearly modulated the local interaction between different spectral components. A comparison of the two experimental conditions revealed a stimulation-related decrease in coupling strength between delta and theta in frontal areas. This coupling seems not to be a target-channel-specific effect of the stimulation on the delta-theta coupling because it overlaps with the area of the strongest default coupling both in STIM and SHAM conditions. Taking into account that the delta and theta power enhancement contains both phase-locked (visible on ERP) and non-phase locked components (visible on TFR) the decrease in coupling between delta and theta in frontal sites might be explained by the altered interaction of these two components of the stimulation-related response. Our subtraction analysis of the evoked and induced components of delta and theta responses on coupling strength revealed that the stimulation-related phase-locked delta component primarily mediates the frontal decrease of coupling. One possible explanation for this observation might be that the slow part of the evoked K-complex propagating frontally is interacting with ongoing theta oscillations. Intracranial recordings could shed light on the interaction between delta and theta in baseline sleep and during stimulus presentation.

In contrast to non-target-electrode-specific delta-theta coupling changes, we found a region-specific effect for coupling between delta and sigma activity. The effect was observed for the area close to the target electrode, with a tendency to spread to the back of the head. This spread might be due to the traveling of slow waves, which primarily propagate from front to back[54]. An explanation of the local nature of the coupling between delta and sigma activity after stimulation might be related to the fact that slow waves are often local and not phase-locked across different brain areas[55]. Given that the auditory stimuli are consistently time-locked to the up-phase in a specific (target) region, the effect might be restricted to that area. Consequently, the stimulation-evoked sigma increase peaked simultaneously with the next up-phase approximately one second after stimulation in that same area.

We also found stimulation-related changes in the preferred phase of sigma oscillations. Following stimulation, the sigma power peak occurred slightly prior to the positive peak of the slow wave, whereas in the SHAM condition power peaked slightly after the positive peak. This is consistent with the observation that external stimulation evokes fast spindles, which are known to appear earlier in the delta cycle compared to slow spindles. Enhancement in the sigma band was also distributed in parietal areas, confirming previous findings on the topography of fast spindles[49,56].In summary, we observed two distinct effects of closed-loop auditory stimulation on cross-frequency coupling for two frequency bands involved in memory consolidation during sleep. In many previous studies, transcranial direct current stimulation was used to demonstrate the relevance of the delta-sigma coupling for the overnight memory formation[24,57,58]. The global sigma boost observed one second after auditory stimulation could provide the optimal time-window for memory replay, with content specificity encoded in the gamma band. In a recent iEEG study, it was shown that the sequences of population firing peaks across widespread cortical regions in the gamma band recorded during wake are re-emerging during NREM sleep. Occurrences of these “motifs” were coupled to down-to-up slow-wave transition, spindles and hippocampal sharp-wave ripples[16]. Along these lines, closed-loop auditory stimulation targeting slow waves in a particular brain region could facilitate the memory replay with particular content by promoting better local coordination between delta, sigma, and nested gamma. In future studies, it will be important to clarify if closed-loop auditory stimulation might improve region-specific memory consolidation.

Study limitations and outlook

One important limitation of our study is the small sample size. However, we believe that the sample size was adequate to demonstrate the impact of closed-loop auditory stimulation on the cross-frequency coupling, as reported effect sizes were within a range from medium to large[59].

In the current study, the stimulation-related change in delta-sigma coupling was topographically restricted. This led us to the assumption that the location of these changes is related to the position of the target electrode. Future studies with different locations of the target electrode would be needed to further clarify the causal association between the location of the electrode used for slow-wave detection and local changes of cross-frequency coupling following the stimulation.

Moreover, it remains to be explored if closed-loop auditory stimulation differently affects slow-wave activity and cross-frequency coupling depending on homeostatic sleep pressure. In the current study, we could not address this question, as stimulation was applied only during the first 2–3 sleep cycles.

In previous studies, it has been demonstrated that at initial stage of the auditory evoked response, it is possible to reconstruct the modality-specific activity of the primary sensory cortices[43,39]. However, our power and cross-frequency coupling analysis necessitated a window of 2 s after stimulus onset and, therefore, modality-specific activity might be masked by the large-amplitude, modality-nonspecific activity originating from associative brain areas.

Conclusion

Our data suggest that closed-loop auditory stimulation might provide a minimally invasive way to modulate the relationships between the main sleep rhythms in a target region while not altering sleep architecture. This opens the way to a large range of future implementations of closed-loop auditory stimulation in different clinical populations. For instance, assessment and local manipulation of cross-frequency coupling during NREM sleep could shed light on the pathophysiological mechanisms underlying sleep-related learning problems observed in clinical populations such as schizophrenia[60,61] or epilepsy[62].

Methods

Participants

Nine right-handed subjects (23 ± 1.3 y.o., 4 females) participated in the study. Participants took notes on their daily activities, sleep and wake time in a sleep diary and were actigraphically monitored (Actiwatch Type AWL from Cambridge Neurotechnology, CamNtech, Cambridge, UK or Geneactiv) for the 7 days (range from 4–8 days) preceding the experiment. Participants did not take any medication at the time of the experiment and were required to refrain from alcohol and caffeine 24 h prior to each experimental session. The sleep schedule in the laboratory was organized in accordance with individual sleep habits. Written informed consent was obtained prior to participation. The study was approved by the local ethics committee (Kantonale Ethikkommission Zürich, KEK-ZH) and performed according to the Declaration of Helsinki.

High-density sleep electroencephalography

Full-night sleep was recorded by using high-density EEG (Electrical Geodesics Sensor Net for long-term monitoring, 128 channels, referenced to a vertex electrode, sampling frequency 500 Hz). Submental electromyographic and electrooculographic data were collected for visual sleep scoring. Two additional electrodes (gold, Grass Technologies, West Warwick, RI, USA) were attached to the earlobes, which served as reference electrodes for the online closed-loop auditory stimulation. After adjusting the net to the vertex and the mastoids, all electrodes were filled with an electrolyte gel to ensure the maintenance of good signal quality throughout the night. EEG electrode impedances were below 50 kΩ. Impedances were below 20 kΩ for the electrodes used for online closed-loop auditory stimulation: submental electrodes, electrodes on the earlobes, and the electrode used for slow-wave detection (C4).

Experimental protocol

Subjects participated in two sessions carried out on different nights separated by 1 week: (i) non-stimulation (SHAM); (ii) closed-loop auditory stimulation in the up-phase (up STIM) of real-time detected slow waves. The order of the STIM and SHAM sessions was randomized and counterbalanced across subjects. In the STIM condition, closed-loop auditory stimulation (STIM) was applied precisely time-locked to the up phase of sleep slow waves detected in the C4 electrode. Real-time slow-wave detection during N2 and N3 was performed as described in Fattinger et al., 2017[23]. The signal of the target electrode C4 was band-pass filtered (Butterworth 0.5–2 Hz) and re-referenced to the mean value of the earlobe electrodes. Auditory stimuli (50 ms bursts of 1/f pink noise, 50 dB with a rising/fall time of 5 ms, inter-stimulus interval ≥ 2 s) were played whenever the EEG signal crossed a default threshold of + 30 μV. The time delay of the stimulus presentation was ∼30 ms. Stimulation was immediately stopped in case of an arousal and if the amplitude of the continuously monitored submental EMG exceeded a given threshold set by the experimenter. Initially, the EMG threshold was defined after 10–15 min of stable N3 sleep in the beginning of the first sleep cycle as an average maximum envelope of the root mean square calculated using a sliding 2 s window. The EMG threshold might have been adapted by the experimenter later during the stimulation period.

Importantly, stimulation was applied only during the first 2–3 sleep cycles, the period with highest levels of SWA. There was no consistent pattern in the distribution of stimuli across the cycles (see supplementary Figure 3).

Data analysis and statistics

EEG data were analyzed using Matlab R2014a, using custom-written scripts and the FieldTrip toolbox (https://fieldtrip.fcdonders.nl/[63]).

Sleep scoring

For sleep scoring, the EEG data were bandpass filtered (FIR filter, 0.5–50 Hz), re-referenced to the earlobes and downsampled to 128 Hz. The EEG was visually scored by two trained individuals for sleep stages Wake, N1, N2, N3 and REM sleep at frontal, central and occipital electrodes (20 s epochs) based on American Academy of Sleep Medicine standard criteria[64].

EEG preprocessing

Data were bandpass filtered (FIR filter, 0.5–45 Hz, Hamming window, mirror padding: 3 s) and preprocessed in overlapping 10 s epochs, with a 6 s prestimulus interval used for baseline correction. Only trials belonging to stages N2 and N3 were considered for the further analyses. In the SHAM condition, where no stimuli were applied, slow-wave up-phase detection was performed offline. Epoched data were downsampled to 128 Hz to speed up processing. Channels located on the earlobes and on the face below the front were excluded from further analysis. After visual inspection, noisy channels and trials containing strong movement artifacts were removed (on average, no more than 5 sensors and 5.8 trials were removed per one recording. Removed channels’ values were interpolated using spherical interpolation. The cleaned data were re-referenced to an average value across all 118 channels. Next, for a consistent phase relationship between the auditory stimuli and the slow waves, we selected only the trials where the stimulus was applied during the rising part of the positive half-wave. For this, we filtered the signal between 0.5–2 Hz and we calculated the instantaneous phase of each slow wave at stimulus onset after applying a Hilbert transformation. The final number of trials included in the analysis did not differ between the two conditions (see supplementary table 2). To compare the number of trials belonging to N2 and N3 epochs, we calculated the percentage of N3 trials considering N2 + N3 as 100%, as only these trials were preselected on the very first step of preprocessing. Again, there was no significant difference between STIM and SHAM conditions.

ERP analysis

Averaged ERP waveforms were computed within each participant and condition. Peak amplitudes were quantified as the average amplitude (± 10 ms) around the local minimum (in case of N550) or maximum (in case of P900) occurring within the timeframe of interest (0.3– 1.3 s for N550 and 0.7–2 s for P900) post-stimulus onset. The time window chosen for the analysis of N550 and P900 was based on visual inspection of the grandaverage waveforms and evidence from previous studies[43].

Frequency analysis

To identify global spectral responses evoked by acoustic stimulation in the time-window of 2 s after stimulus onset, we performed time–frequency analyses separately for total and induced power changes. Total power contains both evoked power (time-locked and phase-locked response to stimulus presentation) and induced power (time-locked but non-phase-locked response). To compute the total power values in STIM and SHAM conditions, single-trial time–frequency representations (TFRs) of power were calculated by FFT with a Hanning taper applied through an adaptive time window of 5 cycles for each frequency between 0.75–40 Hz (ΔT = 3/f for 0.75–3 Hz and ΔT = 4/f for 3–40 Hz) in steps of 0.008 s. Power estimates in the time–frequency domain were averaged across all EEG channels and trials separately for each participant. Through visual inspection of the TFR plot we identified three main spectral contributions in delta (0.75–3 Hz), theta (4–8 Hz) and sigma (13–17 Hz) bands and described their spatiotemporal evolution through topographies. To obtain induced power values, the non-phase-locked activity was calculated by subtracting the ERP from each trial (separately for each condition) and then performing time–frequency analysis on the residual time series.

Cross-frequency coupling

Given the presence of distinct spectral components modulated by the acoustic stimulation, we quantified their relation in terms of cross-frequency coupling strength and preferred coupling phase.

Modulation index (MI)[65] indicating the association between phases in the 0.75–5 Hz (delta) range and power in the 3–40 Hz band was calculated for each 2-s epoch using fieldtrip toolbox[63]. This method is based on the following: under the null-hypothesis of no systematic cross-frequency coupling, the observed amplitude distribution of fast frequency over phase bins of slow frequency does not deviate from the uniform distribution and MI is zero. If the amplitude distribution over phase bins significantly deviates from uniform, the MI will deviate substantially from zero, indicating cross-frequency coupling. For this analysis, we used 18 phase bins (20 degrees each) following[66]. Importantly, coupling strength as assessed by MI depends on the signal-to-noise ratio, which might be confounded by power differences in STIM and SHAM conditions. Therefore, for every subject and channel, we constructed a null distribution of MI by repeatedly (n = 200) shuffling the delta phase timeseries with respect to the fast frequency amplitude timeseries, and recalculating the MI for each iteration (similarly to[5]). This null distribution was used to Z-transform the empirical MI, resulting in a zMI that is independent of power differences. zMI was calculated for all 118 channels included in the power analysis. zMI of 9 representative channels positioned according to the standard 10–20 system were used to demonstrate 2d comodulograms, indicating the cross-frequency coupling in the frequency domain.

The preferred phase analysis was performed in the time-window of 2 s after stimulus onset. After applying a Hanning window (parameters similar to the frequency analysis described above) and obtaining complex Fourier coefficients for the low (delta) and the high (theta or sigma) frequency bands, we computed the absolute amplitude of the high frequency signal in each epoch and sorted the resulting values according to the instantaneous phase of delta oscillation in 180 bins. The preferred phase analysis was performed for 9 representative channels used for comodulograms.

Statistics

To correct for multiple comparisons when examining statistical differences between two conditions, we used a nonparametric clustering procedure[67,68]. First, independent samples t-tests were computed (2-tailed, p < 0.05) for the difference between STIM and SHAM conditions for all tiles in the time–frequency domain. Next, significant neighboring time–frequency tiles were clustered if they showed the same direction of effect. To assess statistical significance of each cluster, a cluster-level test statistic was calculated by computing the sum of all t-values in the cluster. The significance of each cluster-level statistic was estimated by comparing the cluster-level test statistic to a reference permutation distribution derived from the data. The reference distribution was obtained by randomly permuting the data 5,000 times. The cluster p-value was estimated as the proportion of the elements in the reference distribution exceeding the cluster-level test statistic. An analogous procedure was used for the statistical comparison of zMI frequency-frequency plots. Statistical significance of topographical changes was tested in a similar manner: clustering was applied by pooling together neighboring electrodes.

Two-way repeated measures analysis of variance (ANOVA) was used to assess the influence of factors Location (frontal, central, posterior) and Condition (STIM vs. SHAM) on the preferred delta phase separately for theta and sigma bands. A three-ways repeated measures ANOVA with factors Condition (STIM vs. SHAM), EvokedDelta (present or subtracted) and EvokedTheta (present or subtracted) was used to assess the contribution of time-locked components in the coupling of delta and theta. We used the Greenhouse–Geisser correction to estimate the p-values. The level of significance was set to p < 0.05.

The circular Rayleigh test was used to determine whether circular coupling phase distributions deviate from the uniform distribution. All statistical tests were two-sided.

Supplementary information