Co-ordination of brain and heart oscillations during non-rapid eye movement sleep.

Oscillatory activities of the brain and heart show a strong variation across wakefulness and sleep. Separate lines of research indicate that non-rapid eye movement (NREM) sleep is characterised by electroencephalographic slow oscillations (SO), sleep spindles, and phase-amplitude coupling of these oscillations (SO-spindle coupling), as well as an increase in high-frequency heart rate variability (HF-HRV), reflecting enhanced parasympathetic activity. The present study aimed to investigate further the potential coordination between brain and heart oscillations during NREM sleep. Data were derived from one sleep laboratory night with polysomnographic monitoring in 45 healthy participants (22 male, 23 female; mean age 37 years). The associations between the strength (modulation index [MI]) and phase direction of SO-spindle coupling (circular measure) and HF-HRV during NREM sleep were investigated using linear modelling. First, a significant SO-spindle coupling (MI) was observed for all participants during NREM sleep, with spindle peaks preferentially occurring during the SO upstate (phase direction). Second, linear model analyses of NREM sleep showed a significant relationship between the MI and HF-HRV (F = 20.1, r2  = 0.30, p < 0.001) and a tentative circular-linear correlation between phase direction and HF-HRV (F = 3.07, r2  = 0.12, p = 0.056). We demonstrated a co-ordination between SO-spindle phase-amplitude coupling and HF-HRV during NREM sleep, presumably related to parallel central nervous and peripheral vegetative arousal systems regulation. Further investigating the fine-graded co-ordination of brain and heart oscillations might improve our understanding of the links between sleep and cardiovascular health.

INTRODUCTION

Numerous studies corroborate important functions of sleep for the brain, such as synaptic scaling (Tononi & Cirelli, 2014), metabolic clearance (Ding et al., 2016), and overall mental health (Baglioni et al., 2016). Other studies address the relevance of sleep for peripheral organs, such as cardiovascular regulation (Kuo, Chen, Hsu, & Yang, 2016; Penzel et al., 2016) and diseases (Lane et al., 2019; Pawlowski et al., 2017), but the complex interplay of central nervous system (CNS) and peripheral functions need to be further specified.

The transition from wake to sleep is characterised by deactivation of the main arousal system in mammals, the ascending reticular arousal system (ARAS), and associated changes in brain oscillatory activity. A major characteristic of light non‐rapid eye movement (NREM) sleep (N2) are thalamo‐cortical oscillations in the form of sleep spindles (sigma frequency range, 12–16 Hz). During deep NREM sleep (slow‐wave sleep [SWS], N3), increasing synchronisation of neural activity is reflected by the occurrence of cortical slow oscillations (SO; 0.5–2 Hz) and high amplitude delta frequency waves (1–4 Hz; Andrillon et al., 2011; Molle, Marshall, Gais, & Born, 2002). SOs have been observed to orchestrate other oscillations, such as sleep spindles (Contreras & Steriade, 1995). This co‐ordination can be quantified by the strength of phase–amplitude coupling (PAC) between these oscillations (SO–spindle phase–amplitude coupling). More specifically, PAC describes the temporal interaction between the phase of a low frequency rhythm (i.e. SOs) and the amplitude of a higher frequency component (i.e. sleep spindles; Canolty et al., 2006; Helfrich, Huang, Wilson, & Knight, 2017; Voytek et al., 2010). This means that the power envelope (amplitude) of the spindle’s changes with the phase of the SOs.

Recent work suggests that PAC between SOs and sleep spindles promotes brain functions, such as attention, decision‐making(Canolty & Knight, 2010) or the consolidation of memory traces (Helfrich, Mander, Jagust, Knight, & Walker, 2018; Mikutta et al., 2019; Staresina et al., 2015) and underlying neural plasticity (Maier et al., 2019). Of note, not only the strength of SO–spindle coupling, but also the temporal orchestration (i.e. the position of the amplitude during the phase) appears to be critical for information processing during sleep (Helfrich et al., 2018, 2019; Winer et al., 2019). For example, the position of the spindle peak relative to the SO phase (i.e. phase direction) determined overnight memory consolidation, with individuals showing the spindle peak on the up‐state of the SO having the best outcome (Helfrich et al., 2018). Moreover, the phase direction was altered (spindle peak before SO up‐state) in individuals with diminished frontal lobe volumes (Helfrich et al., 2018).

Until now, PAC has primarily been related to brain functions. Yet, in a broader perspective, enhanced PAC might reflect deactivation of arousal networks and reduced susceptibility to external stimuli during sleep (Dang‐Vu, McKinney, Buxton, Solet, & Ellenbogen, 2010). We reasoned that such a deactivation of CNS activity might go along with deactivation of peripheral body functions, as indexed by a shift in autonomic nervous system (ANS) function from sympathetic towards parasympathetic vagal activity. Sympathetic activity relates to the fight‐flight or freeze responses. Parasympathetic activity, in turn, reflects the component of the ANS that promotes rest‐related activities, such as “rest and digest” or “feed and breed” activities (Wehrwein, Orer, & Barman, 2016). Prior work identified heart rate variability (HRV) as a reliable index of parasympathetic activity (Jung, Lee, Jeong, & Park, 2017; Camm et al., 1996). HRV embraces the concept that the time between adjacent heartbeats differs in length. This variability is often quantified by computing power in spectral frequency bands. The role of the low‐frequency HRV (LF‐HRV: 0.04–0.15 Hz) remains unclear due to the underlying mixed ANS influences (Goldstein, Bentho, Park, & Sharabi, 2011; Martelli, Silvani, McAllen, May, & Ramchandra, 2014; Rahman, Pechnik, Gross, Sewell, & Goldstein, 2011). In contrast, the physiological mechanisms underlying rapid changes in HRV are well characterised. High‐frequency HRV (HF‐HRV: 0.15–0.40 Hz) reflects cardiac parasympathetic activity (Chess, Tam, & Calaresu, 1975; Pomeranz et al., 1985; Thayer, Hansen, Saus‐Rose, & Johnsen, 2009). Along these lines, greater parasympathetic activity, indexed by HF‐HRV, is observed during NREM sleep compared to wakefulness (Carrington et al., 2005; Grimaldi et al., 2019; Mendez, Bianchi, & Cerutti, 2004; Nazeran et al., 2006).

The present study aimed to assess, to our knowledge for the first time, the interaction of brain (PAC) and heart (HRV) oscillations during sleep. To this end, we conducted novel analyses of PAC and its relation to HRV in an existing dataset of sleep laboratory nights in healthy participants. Particularly, we tested the hypothesis that brain and heart oscillations during NREM sleep would be co‐ordinated, operationalised as a positive correlation between the strength (modulation index [MI]) of SO–spindle PAC (phase: SO, 0.5–2 Hz; amplitude: sigma activity, 12–16 Hz) and HF‐HRV, and a circular correlation between the phase direction of SO–spindle phase–amplitude coupling and HF‐HRV. In a post hoc exploratory analysis we estimated the interaction between PAC and MI as described above during the first three sleep cycles.

METHODS

Participants

Data from 45 healthy participants (22 males, 23 females, aged 21–60 years, mean [SD] age 37.2 [12.3] years) were analysed. The number of participants was based on a power analysis (G*Power) with a significance level alpha of 0.05, a power of 0.8 and an effect size f 2 of 0.3 for linear regression analyses. Participants were consecutively drawn from existing datasets from the sleep laboratory of the Department of Psychiatry and Psychotherapy of the Medical Center, University of Freiburg, Germany (Feige et al., 2008; Holz et al., 2012). None of the participants had a history of sleep, neurological, psychiatric, or other somatic disorders according to a clinical interview by an experienced physician. Participants with an apnea–hypopnea index of >5 events/hr were excluded. All participants were right‐handed, non‐smokers, and free of any medication or drug use. All participants provided written informed consent. The work was approved by the Ethics Committee of the University of Freiburg, Germany, and was carried out in accordance with The Code of Ethics of the World Medical Association (Declaration of Helsinki).

Polysomnographic recordings

As described in detail previously (Holz et al., 2012), polysomnography (PSG) was recorded from 11:00 p.m. to 7:00 a.m. in the sleep laboratory of the Department of Psychiatry and Psychotherapy of the Medical Center, University of Freiburg, Germany. Data were drawn from one baseline night after one adaptation and screening night (including respiratory signals and oxygen saturation). Sleep was recorded using 24‐channel Sagura electroencephalography (EEG)‐PSG systems (Sagura Polysomnograph 2000, Sagura Medizintechnik GmbH). All signals were recorded synchronously with the same system and stored in a single multi‐channel European Data Format (EDF) file. Sleep recordings were scored by experienced raters based on standard criteria (Rechtschaffen & Kales, 1968.). The setup included the EEG electrodes C3–A2 (sampling rate 256 Hz), submental electromyogram, vertical and horizontal electro‐oculogram and electrocardiogram (ECG, Lead II Einthoven configuration, sampling rate 512 Hz). For details on pre‐processing of the EEG data please refer to Holz et al. (Holz et al., 2012). All further analysis steps used the same artefact‐free overlapping segments of EEG and ECG data. The following PSG parameters were calculated: sleep latency, total sleep time, rapid eye movement (REM) sleep latency, SWS latency, and, with regard to sleep architecture, percentages of the different sleep stages referred to total sleep time.

Selection criteria for EEG and ECG data segments

Each participant had to have >120 min of overall artefact‐free EEG data during NREM sleep and the equivalent ECG data set of each participant had to have <5% artefacts using the automatic artefact correction algorithm implemented by Kubios 3.4 (premium version, https://www.kubios.com). In order to maximise the quality and amount of data included into further analyses, 1‐min data segments of the ECG and EEG data had to be without sleep stage transitions and without artefacts in order to ensure that EEG and HRV analyses (see chapter Heart rate variability analysis) could be performed on the same segments. The choice of 1‐min segments length is a trade‐off between the two signal components (ECG and EEG).

Phase–amplitude coupling analysis

All analyses were done using Matlab 2019 (https://www.mathworks.com). In a first step, to explore the phase–amplitude coupling between low‐frequency phases and a high‐frequency amplitude in our data, the comodulogram averaged across all NREM sleep stages and across all participants was calculated. The comodulogram is a data‐driven approach indicating the strength of coupling between the phase of a low‐frequency part of the signal and the amplitude of a high‐frequency part, with the x‐axis representing the low and the y‐axis the amplitude‐modulated high frequency (Figure 1d and additional Figure 1b). For the comodulogram the MI (Tort, Komorowski, Eichenbaum, & Kopell, 2010) was estimated, as described previously, (Maier et al., 2019; Mikutta et al., 2019) on 1‐min data segments. Data segments were chosen from artefact‐free (both EEG and ECG) epochs without sleep stage transition to ensure that both EEG and HRV analyses (see chapter Heart rate variability analysis) could be performed on the same segments. We used finite‐duration impulse response filters to band‐pass filter the data (forward/backward, filtfilt, Matlab, signal processing toolbox, https://www.mathworks.com/products/signal.html). Finite impulse response (FIR) filters used an order equal to three cycles of the lower cut‐off frequency. In brief, the MI is defined as an adaptation of the Kullback–Leibler distance (a function that is used to infer the distance between two distributions) and calculates how much an empirical amplitude distribution‐like function over phase bins deviates from the uniform distribution. In detail, the MI was computed as follows: The phase of the SO was divided into B = 18 bins (each covering 20 degrees). The mean power of the spindle frequency oscillations envelope within each bin was estimated. By combining the individual bins, a phase–amplitude histogram (PAH) was computed. We compared the PAH with a uniform distribution (UD) using the Kullback‐Leibler distance (DKL):

FIGURE 1

(a; top) Representative raw ECG data from one participant. Red asterisks denote the peak of the QRS complex, R wave. The raw data show that the R–R interval changes. (bottom) Time between R–R intervals for the ECG shown in A (Top). (b; left, top and bottom) Frequency domain analysis: Power spectral density (PSD) of the autoregression analysis of the R–R intervals. The default values for the bands are VLF: 0–0.04 Hz (black), LF: 0.04–0.15 Hz (blue), and HF: 0.15–0.4 Hz (red). The example depicts data with high HF frequency power. (c; top) An example of the EEG raw signal showing a slow wave event at approximately 2 s. (c; bottom) The same EEG data filtered in the SO range (0.5–2 Hz, red) and spindle‐range (12–16 Hz, blue). The highest amplitudes of the spindle occur in the up‐state of the SO. (d) Comodulogram between lower frequency (0.5–6 Hz, 0.5 Hz steps) phases and faster frequency (8–35 Hz, 1‐Hz steps) amplitudes during NREM sleep. (e) Preferred phase of the 45 participants on a circular plot. 0/pi (0/360°) equals the up‐state of the SO wave. The red line depicts the mean vector over all participants. ECG, electrocardiogram; EEG, electroencephalography; HF, high frequency; LF, low frequency; NREM, non‐rapid eye movement; SO, slow oscillation; VLF, very low frequency

The DKL was normalised to calculate the MI:

For the comodulogram, the MI between low frequencies of 0.5–6 Hz in 0.5‐Hz steps and faster frequencies of 8–35 Hz in 1‐Hz steps was calculated. As for this exploratory analysis no specific events, but frequency ranges were defined, we refer to the standard frequency range specifications (that is, SO and sigma activity).

In a second step, the individual strength and phase direction of SO–spindle phase–amplitude coupling during NREM sleep of each participant was calculated, as described previously (Mikutta et al., 2019). In detail, for this SO event‐locked approach, we detected the SO events using the algorithm described in Esser et al. (Esser, Hill, & Tononi, 2007) and cut out SO‐locked events using 1.5 s before and after each SO peak. We z‐normalised events per participant prior to all analyses. Mean and standard deviations (SDs) were computed using the unfiltered event‐locked SO in every participant to avoid spurious coupling (Cole & Voytek, 2017). Furthermore, we defined spindle events by extracting the analytical amplitude after applying a Hilbert transform. We smoothed the amplitude with a 200 ms moving average. Then the amplitude threshold was set at the 75% percentile (amplitude criterion) and only events that exceeded the threshold for 0.5–3 s (time criterion) were accepted. For further analyses, we selected only SO events, with a spindle within a ±1.5 s window around the SO peak. Due to this event selection we refer to SO and spindle events in this analysis. For estimating the individual strength of SO–spindle coupling, the MI for each participant and event, data was band‐pass filtered within the SO (0.2–5 Hz) and spindle (12–16 Hz) frequency bands (Figure 1a,b). For avoiding filter edge artefacts, we discarded the first and last 0.5 s of each SO epoch. Computation of the MI again followed the procedure outlined above. The individual MI for each participant is reported as average over all trials. For estimating the individual phase direction of SO–spindle coupling, the instantaneous phase of the SO band (0.5–2 Hz) for all events was extracted. After that, the instantaneous amplitude of the same events within the sigma band (12–16 Hz) was extracted. Next, we detected the maximal peak of the sigma amplitude and the corresponding phase angle of the SO in every subject and epoch. Then, the preferred phase was computed by calculating the mean sigma peak relative SO phase angle of all SO epochs within each participant using the circ‐mean function of the CircStat toolbox (15). For linear modelling, we converted this circular data (preferred phase) into the cartesian space. In order to keep out the vector length (strength of interaction) and solely focus on the phase angle, we used one instead of the vector length as multiplicator for getting the Real and Imaginary part in the artesian space: Real = sin(preferred phase) *1; Imaginary = cos(preferred phase) *1.

Heart rate variability analysis

All analyses comprising HRV parameters represent novel, previously unpublished analyses. The ECG data were acquired at a 512‐Hz sampling rate using a Lead II Einthoven configuration. Data were down sampled to 500 Hz using Matlab 2019a. Analysis of HRV was conducted according to the Task Force of the European Society of Cardiology guidelines (Camm et al., 1996). However instead of 5‐min excerpts, we only used 1‐min data excerpts as described above. The 1‐min windows were chosen in order to maximise the amount of artefact‐free EEG and ECG data for the analyses. Recent research suggests that 1‐min time windows are sufficient for providing reliable HRV results (Koenig et al., 2017).

The HRV data was characterised by providing absolute measures of total power, HF power in ms2 units and the heart rate (HR) in beats per min. Furthermore, HF‐peak frequency as a surrogate marker for the breathing rate is provided. Several studies confirmed the reliability of ultra‐short HRV frequency analysis (Lischke et al., 2018; Salahuddin, Cho, Jeong, & Kim, 2007; Seps, Beckers, & Aubert, 2002). Time windows were averaged across all NREM sleep stages, which is preferable for ultra‐short HRV recordings (Schroeder et al., 2004). HRV of the R‐wave series derived from the ECG were analysed using MATLAB based Kubios HRV Analysis Software 2.0 (Tarvainen, Niskanen, Lipponen, Ranta‐Aho, & Karjalainen, 2014). An autoregressive algorithm (smoothing prior Lambda 500, artefact correction threshold low) was used to compute the spectral power in the HF (0.15–0.40 Hz) band, in ms². All values are reported as means over all 1‐min artefact‐free NREM epochs of each participant (as defined in Selection criteria for EEG and ECG data segments).

Sleep cycle analysis of SO‐spindle coupling and HF‐HRV

As both variables, that is SO–spindle coupling and HF‐HRV, are linked to time windows, we wanted to assure that the observed effects were stable across different time scales by using the different sleep cycles as a confounding factor within linear modelling. We explored changes of HF‐HRV, MI and the relationship thereof during the night, we separately analysed NREM sleep data of the first three sleep cycles. Participants included in the sleep cycle analyses had to have >20 min of artefact‐free data in each sleep cycle.

Statistical analysis

Descriptive values for the PSG are provided as means (SDs). Furthermore, we provide descriptive values for the comodulogram, individual MI and HF‐HRV as means and SDs, if not stated otherwise. For testing the individual strength (MI) of SO–spindle coupling for significance, we performed a surrogate control analysis by shuffling the spindle amplitudes and SO phases over trials and thereby creating 200 surrogate MI values for estimating a significance threshold (p < 0.05) for each participant. Before we computed the preferred phase, we controlled for a non‐uniform distribution of phase directions of SO–spindle coupling within each participant, using a Rayleigh test for non‐uniformity of circular data. To test the primary hypothesis of a linear relationship between the strength of SO–spindle coupling and HF‐HRV during NREM sleep, we fitted a linear model with MI as predictor and HF‐HRV as outcome parameter. An F test of overall significance (comparison to an intercept only model) was conducted and the explained variance (R 2) and the B value (beta; i.e. coefficients or differences between factor levels) of the fitted model are provided. To test the second primary hypothesis, we used circular‐linear correlation to estimate the correlation between HF‐HRV (linear) and the preferred phase (circular). The r value and corresponding p value is reported.

To control for the stability of our results for different time windows, we computed individual MI, phase direction and HF‐HRV values for each of the first three sleep cycles and compared those using one‐factorial analyses of variance (ANOVAs) with the within‐subject factor Sleep cycle for MI and the cartesian co‐ordinates (real and imaginary part) and HRV‐HF. For testing for differences in the preferred phase, we used a Watson–Williams multi‐sample test for equal means. For post hoc comparisons, paired Student’s t tests were conducted. The p values of the paired Student’s t tests were corrected for multiple testing using false discovery rate (FDR) correction (Benjamini and Hochberg method). To test for the impact of the factor sleep cycle, we fitted a repeated measure linear model using the individual MI value of each of the three sleep cycles (three values per participant) as input parameters and the HF‐HRV values with cycle as confounding variable. We report the results as F test of overall significance. For univariate effect, we report F and p values, the explained variance (r 2) of the predictor variables as well as B values (betas) of the linear model, i.e. coefficients or differences between factor levels. The p values of within‐subject parameters were Greenhouse–Geisser corrected. Furthermore, we fitted a repeated measure linear model using the individual real and imaginary part cartesian co‐ordinate values of each of the three sleep cycles (three values per participant) as input parameters and the HF‐HRV values with cycle as confounding variable. We report the results as outlined above. In order to control for alternative explanations, we tested for confounding factors (duration of NREM overall, SWS, number of SO events, EEG delta power, EEG sigma power, age) and concerning HRV (HF peak frequency [as a proxy for breathing frequency], SD of the starting times of 1‐min segments over the night [as proxy for the overnight distribution of the segments], number of R‐peaks) by increasing the linear models with these confounding parameters as predictors.

RESULTS

Polysomnography

For the primary analyses of NREM sleep of the entire night, we included all 45 participants. For the exploratory sleep cycle analyses, we included 36 participants meeting the inclusion criteria (>20 min of artefact‐free ECG and EEG data in each sleep cycle). In brief, PSG recordings demonstrated a mean (SD) sleep latency of 31.3 (20.1) min, total sleep time of 470.6 (23.3) min, NREM sleep of 271.4 (65.4) min, SWS of 36.2 (26.7) min, REM sleep latency of 91.1 (25.7) min, REM sleep of 85.1 (28.5) min and, with regard to sleep architecture, 67% of NREM sleep (9% of SWS) and 18% of REM sleep referred to the total sleep time. For the current analyses, a mean (SD) 188.0 (53.9) min of NREM sleep including 28.3 (22.4) min of SWS were used. Data for single sleep cycle is provided in Table S1.

Strength of phase‐amplitude coupling

For the calculated comodulogram of the whole night’s NREM sleep episodes, the highest strength of coupling was found between the standard SO phase (0.5–2 Hz) and the sigma (12–16 Hz) amplitude (Figure 1d), as previously reported in our own work (Mikutta et al., 2019) and that of other groups (Helfrich et al., 2018).

Second, in order to estimate the individual strength of the coupling between SO and spindle activity in each participant, we identified a mean (SD) of 1837 (725) SO events with a sleep spindle in ±1‐s time‐window around the peak during NREM overall sleep. Replicating prior observations (Helfrich et al., 2018; Klimesch, 2013), a strong coupling for the standard SO phase (0.5–2 Hz) and spindle frequency range amplitude (12–16 Hz) was observed (NREM sleep: mean [SD] MI across all participants 0.0119 [ 0.0059]). Individual values are presented in Figure 2 and additional Figure 2. For all participants individual MI values for this SO phase and spindle frequency range were calculated using the surrogate MI values and yielded to be significantly higher (p < 0.05) than the threshold based on compared surrogate data, indicating a significant phase–amplitude coupling between the SO and spindle frequency range in each participant.

FIGURE 2

(a) HF‐HRV (ms2) as a function of MI. We found a significant correlation (r = 0.56, p = 0.001), indicating that a stronger interaction between spindle and SO goes along with a higher parasympathetic activity. HF‐HRV (ms2) as a circular‐linear function of the mean phase direction during NREM sleep. (b) The tendency towards a significant circular‐linear correlation between HRV‐HF and the phase direction (r = 0.35, p = 0.056) indicates that a precise orchestration of spindle peaks and SO up‐states coincides with a higher HF‐HRV values indicating increased parasympathetic activity. The black line indicates a quadratic fit for visualising the circular‐linear correlation. (c) Linear model (adjusted values) between HF‐HRV and MI. The model reveals a significant interaction between HF‐HRV and MI (F [1,42] = 20.1, r 2 = 0.30, B‐value = 5.62, p < 0.001). HF‐HRV, high‐frequency heart rate variability; MI, modulation index; NREM, non‐rapid eye movement; SO, slow oscillation

Phase direction of phase–amplitude coupling

Next, in order to explore the specific temporal organisation between SO and spindle events, we identified the mean coupling direction, that is the SO phase‐angle position of the maximal spindle amplitude within each participant. The Rayleigh test showed a significant non‐uniform circular distribution in every participant (p < 0.01, for all participants). We found 64% of the participants (29/45 participants) having the preferred phase on the up‐state (window ±30° window). Particularly, the highest peaks of the spindle amplitude were located on average at a mean (SD) of −0.010 (0.62) radians during NREM overall sleep, that is during the up‐states.

Heart rate variability

The mean (SD) heart rate during NREM sleep was 50.3 (15) beats/min. The mean (SD) total power was 4921 (1295) ms2 and the HF‐HRV was 2728.3 (698.1) ms2. The mean (SD) HF‐peak frequency was 0.16 (0.08) Hz. Analyses were based on a mean (SD) of 10911 (8263) beats during NREM sleep. Individual HF‐HRV can be derived from Figure 2 and additional Figure 2.

Linear relationship between strength of SO‐coupling and HF‐HRV

The linear regression model between the MI and HF‐HRV during NREM sleep was significant (F [1,42] = 20.1, r 2 = 0.30, B‐value = 5.62, p < 0.001; Figure 2c), indicating a significant linear relationship between SO‐coupling and HF‐HRV (Figure 2a and c).

Circular‐linear relationship between phase direction of SO‐coupling and HF‐HRV

The circular‐linear between preferred phase and HF‐HRV demonstrated a tendency towards significance without attaining statistical significance (r = 0.34, p = 0.056).

Post hoc sleep cycle analyses

In the comodulogram for NREM sleep in each sleep cycle, we observed a decreasing strength of SO–spindle coupling (additional Figure 1b, left panel). For estimating the individual strength of the coupling between SO and spindle activity in each participant for each sleep cycle, we identified a mean (SD) ncycle 1 = 496 (262), ncycle 2 = 485 (211), ncycle 3 = 448 (182) SO events with a sleep spindle in ±1‐s time‐window around the peak. Across sleep cycles, the MI decreased significantly in the third sleep cycle (one‐factorial ANOVA: F [2,105] = 3.78, p = 0.02). For the post hoc Student’s t test the p value was adjusted to 0.02 using the FDR (cycle 1 versus cycle 3: t = 3.60, p = 0.001; cycle 2 versus cycle 3: t = 2.67, p = 0.011; Figure 2a, left panel).

During the first three sleep cycles, the mean (SD) phase direction of SO–spindle coupling was −0.38 (0.14), −0.41 (0.21), and −0.21 (0.09) radians, respectively. There was no significant difference of phase direction between the three sleep cycles (one‐factorial ANOVA: F‐real[2,105] = 0.19, p = 0.82; one‐factorial ANOVA: F‐imagery[2,105] = 0.14, p = 0.86).

The HRV parameters for all three sleep cycles are listed in Table 1. There was no significant difference in HF‐HRV between the three sleep cycles (one‐factorial ANOVA: F [2,105] = 0.16, p = 0.85; additional Figure s1a, right panel).

TABLE 1

Heart rate variability parameters for the first three sleep cycles

	Sleep cycle 1	Sleep cycle 2	Sleep cycle 3
Heart rate, beats/min, mean (SD)	52 (15)	51 (12)	52 (14)
Total power, ms2, mean (SD)	3995 (1376)	4027 (1391)	3857 (1150)
HF power, ms2, mean (SD)	2735 (950.9)	2673 (921.4)	2715 (871.2)
HF‐peak frequency, Hz, mean (SD)	0.13 (0.07)	0.15 (0.09)	0.12 (0.03)
Beats analysed, n, mean (SD)	2850 (2963)	2939 (3714)	2502 (2935)

HF, high frequency.

The preferred phase did no change significantly between the sleep cycles (Watson–Williamson test: F [2,105] = 0.29, p = 0.75; Figure s1c).

Further, for testing the stability of the observed significant relationship between MI over the whole night throughout NREM sleep, the fitted repeated measures linear model with the predictors MI and sleep cycle as within subjects variables and the outcome parameter HF‐HRV was significant (F [1,33] = 19.7, r 2 = 0.23, B‐value = 2.04, p < 0.001; Figure s2a and c). Furthermore, the number of the sleep cycle was not a significant within subject factor (F [1,68] = 1,9, p GG =0.14; Figure s2c), indicating the stability of the relationship between SO–spindle coupling and HF‐HRV for the first three sleep cycles.

The repeated measure linear regression model with the outcome parameter HF‐HRV and the predictors real and imaginary part using the variable cycle as within subject parameter did not show a significant fit (F [2,31] = 1.08, r 2 = 0.06, B‐value = −0.18, p = 0.35), real part (F real = 0.78, r 2 = 0.04, B‐value = −0.02, p = 0.46), and imaginary part (F imag = 5.3, r 2 = 0.25, B‐value = 0.03, p = 0.01; Figure s2b and d). Cycle was not a significant within subject factor (F real = 0.9, p GG = 0.4; F imag = 1.04, p GG = 0.364).

Control analyses

Results of the control analyses can be found in the Supplementary Material.

DISCUSSION

The present study shows, for the first time, to our knowledge, that the strength and precise temporal co‐ordination (phase direction) of sleep slow oscillatory and spindle activity, quantified as PAC, during NREM sleep in the brain goes along with an increased HF‐HRV, an index of enhanced peripheral parasympathetic activity. This finding contributes to and extends the emerging knowledge on orchestrated brain and body functions during sleep.

Previously, the concept of PAC was primarily associated with cognitive functions, e.g. memory consolidation (Canolty & Knight, 2010; Helfrich et al., 2018). Our present results are consistent with recent findings indicating that peripheral cardiac functions (ECG) and EEG patterns show distinct synchronisation patterns (Klimesch, 2013). In detail, it was shown that ECG R‐peaks coincided with up‐states of SO, as well as spindle peaks during SWS (Lechinger, Heib, Gruber, Schabus, & Klimesch, 2015). Our present work extends prior work by demonstrating that the strength of the SO–spindle interaction (MI) correlates with a measure of parasympathetic activity. Therefore, our present finding of strength of coupling and phase direction correlating with autonomic activity might establish a further link between central brain functions and HRV, indicating that a precise temporal orchestration of SO and spindle events coincides with a higher parasympathetic activity. Our present study was not designed to control for directionality of this process. However, recent work assumes a bi‐directionality of communication between neural oscillations and ECG R‐peaks (Grimaldi et al., 2019).

Although it is too early for mechanistic conclusions, we speculate that the strength of SO–spindle coupling might reflect an increased general synchronisation during sleep, which indicates a decreased impact of external stimuli on the brain. This, in turn, might cause an increased parasympathetic activity. This view is supported by studies that observed cortical and peripheral autonomic activity to be macroscopically associated (Jurysta et al., 2003). A body of research indicates that the temporal interaction of cardiac events and delta and sigma power in the EEG predicts memory performance (Naji, Krishnan, McDevitt, Bazhenov, & Mednick, 2019a, 2019b), indicating the general importance of brain (EEG) and autonomous (ECG) activity interaction (Whitehurst, Chen, Naji, & Mednick, 2020).

Specifically, prior work suggested an interaction between the absolute power of most conventional EEG frequency bands and HF‐HRV as a marker for parasympathetic activity (Otzenberger et al., 1998; Otzenberger, Simon, Gronfier, & Brandenberger, 1997). A specific positive association was found between delta power and HRV‐HF.4,49,55 Further studies indicate a relevance of the phase of EEG delta activity, as delta phase shifts precede changes in HF‐HRV activity (Jurysta et al., 2003). Non‐linear HRV analyses methods revealed that the occurrence of SO/delta oscillations interact with non‐linear HRV parameters (Yeh et al., 2013). However, it remains to be further clarified how the rather phasic measure of MI relates to the rather tonic HF‐HRV.

The general physiological bases of PAC are still under debate. It was shown that a disturbed structural integrity of the prefrontal cortex due to ageing is associated with a dissociation of SO and spindles (Helfrich et al., 2018). The prefrontal cortex was as well identified as a key structure for HRV regulation (Thayer et al., 2009). Therefore, we speculate that a physiological origin for the correlation between PAC and HRV might be located in the prefrontal cortex.

The control analyses confirmed the exclusive relationship between HF‐HRV and MI, as none of the conventional (absolute and relative) frequency bands showed a significant correlation. Further, other HRV parameters (SD of normal‐to‐normal R–R intervals [SDNN], root mean square of successive R–R interval differences [RMSSD]) did not show a significant correlation with MI or preferred phase. The additional analysis of SWS did not show any difference to the analysis of NREM sleep and therefore did not change the main outcome.

Although the absolute values of the MI declined during the third sleep cycle, the correlation between MI and HF‐HRV remained significant and was not different to the first and second cycle. We speculate that the relative changes and not the absolute strength of the MI are relevant for the interaction with HF‐HRV.

Major limitations of the present study comprise the rather small but according to power analysis statistical reasonable number of participants and the fact that we re‐analysed the data from an existing dataset. The compromises of analysing methods due to the data structure (EEG and ECG data) limited on the one hand the possible time resolution (low validity of the MI with a low number of trials, 1‐min time‐window for HF‐HRV calculation) and on the other hand prerequisites the need to discard data in order to get the optimal data quality for EEG and ECG data in the respective 1‐min epochs. A further limitation is the selection of HRV variables. Due to the given requirements of the data, we opted for frequency domain HRV analysis. We chose HF ms2 based on the recent literature, indicating HF ms2 to be a reliable parameter of parasympathic activity (Grimaldi et al., 2019). Although we controlled our present results using time‐domain variables (RMSSD), there might be a restriction of these results due to the short analysis window. Another limitation is the analysis of a single electrode. However, recent results indicate that left fronto‐parietal areas covered by C3 show strong coupling activity during sleep (12). Further frontal cortex areas have often been associated with central control of HRV (30).

Further developing the present line of research might have relevant clinical implications. For instance, a lowered HF‐HRV is associated with an increased risk of cardiovascular incidents (Hillebrand et al., 2013). Sleep, and in particular orchestrated brain and heart oscillations, might provide a window for new interventions, such as through non‐invasive brain stimulation techniques, including auditory stimulation (Bonnet & Ardand, 2007; Fehér, 2020; Peppard, Young, Palta, & Skatrud, 2000). In particular, the quality and the fine‐graded orchestration of physiological processes during sleep might be more important than pure sleep duration.

The results of our present study might have relevant clinical implications. Our present findings support the notion that disturbed orchestration of brain oscillations, indicative for enhanced arousal during sleep, has a potential link to cardiovascular diseases (Bonnet & Ardand, 2007; Peppard et al., 2000).

Together, the present study demonstrates an interaction between MI, phase direction and an index of HRV, with higher MI values going along with indices of increased parasympathetic activity. Further analyses of the fine‐graded orchestration of brain oscillations during sleep might have the potential to improve our understanding of the links between sleep and peripheral functions, such as cardiovascular activity and disease. Future studies are needed to explore possible therapeutic agents and interventions to alter PAC and related functions.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Conceptualisation: CM, MW, and CN. Methodology: CM, BF, AA, MW, KS, and JK. Software: CM, BF, AA, and JK. Validation: CN, DR, AA, and KS. Formal analysis: CM, BF, AA, MW, JM, and JK. Investigation: BF and CN. Resources: CN and DR. Writing – Original Draft: CN, CM, MW, EH, CS, and JM. Visualisation: CM, MW, BF, and KF. Supervision: CN and BF.

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Supplementary Material

Click here for additional data file.