Front. Physiol. Frontiers in Physiology Front. Physiol. 1664-042X Frontiers Media S.A. 1511998 10.3389/fphys.2024.1511998 Physiology Original Research Detection of respiratory frequency rhythm in human alpha phase shifts: topographic distributions in wake and drowsy states Kalauzi et al. 10.3389/fphys.2024.1511998 Kalauzi Aleksandar 1 Matić Zoran 2 Suljovrujić Edin 2 Bojić Tijana 2 * 1 Department for Life Sciences, Institute for Multidisciplinary Research, University of Belgrade, Belgrade, Serbia 2 Laboratory for Radiation Chemistry and Physics-030, Institute for Nuclear Sciences Vinča-National Institute of the Republic of Serbia, University of Belgrade, Belgrade, Serbia

Edited by: Alberto Porta, University of Milan, Italy

Reviewed by: Beatrice Cairo, University of Milan, Italy

Raphael Martins de Abreu, Lunex University, Luxembourg

 *Correspondence: Tijana Bojić, tijanabojic@vin.bg.ac.rs 06 01 2025 2024 15 1511998 18 10 2024 09 12 2024 Copyright © 2025 Kalauzi, Matić, Suljovrujić and Bojić. 2025 Kalauzi, Matić, Suljovrujić and Bojić

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

 Introduction

The relationship between brain activity and respiration is recently attracting increasing attention, despite being studied for a long time. Respiratory modulation was evidenced in both single-cell activity and field potentials. Among EEG and intracranial measurements, the effect of respiration was prevailingly studied on amplitude/power in all frequency bands.

 Methods

Since phases of EEG oscillations received less attention, we applied our previously published carrier frequency (CF) mathematical model of human alpha oscillations on a group of 10 young healthy participants in wake and drowsy states, using a 14-channel average reference montage. Since our approach allows for a more precise calculation of CF phase shifts (CFPS) than any individual Fourier component, by using a 2-s moving Fourier window, we validated the new method and studied, for the first time, temporal waveforms CFPS(t) and their oscillatory content through FFT (CFPS(t)).

 Results

Although not appearing equally in all channel pairs and every subject, a clear peak in the respiratory frequency region, 0.21–0.26 Hz, was observed (max at 0.22 Hz). When five channel pairs with the most prominent group averaged amplitudes at 0.22 Hz were plotted in both states, topographic distributions changed significantly—from longitudinal, connecting frontal and posterior channels in the wake state to topographically split two separate regions—frontal and posterior in the drowsy state. In addition, in the drowsy state, 0.22-Hz amplitudes decreased for all pairs, while statistically significant reduction was obtained for 20/91 (22%) pairs.

 Discussion

These results potentially evidence, for the first time, the respiratory frequency modulation of alpha phase shifts, as well as the significant impact of wakeful consciousness on the observed oscillations.

 respiration electroencephalography alpha activity phase shift wake and drowsy mind–body interaction section-at-acceptance Autonomic Neuroscience 1 Introduction

Breathing, one of the fundamental life rhythms, occurs at a frequency of approximately 0.25 Hz in humans at rest and consists of active inspiration and mostly passive expiration (Fleming et al., 2011). For many years, the relationship between brain activity and respiration has attracted researchers’ attention, and different quantitative measures, referring to this relationship, have been applied (Harper et al., 1998; Heck et al., 2017; Herrero et al., 2018; Del Negro et al., 2018; McKay et al., 2003; Nakamura et al., 2023). Due to its nature, respiration could be easily tracked and recorded, while the brain acts as an incomparably richer source of information, resulting in numerous methods applied to monitor its respiration-related activity. In addition to electroencephalography (EEG), most common methods include (BOLD) fMRI (Damoiseaux et al., 2006; Jung and Kang, 2021), intracortical recording (Viczko et al., 2014; Väyrynen et al., 2023), or magnetoencephalography (MEG) (Hsu et al., 2020; Kluger and Gross, 2020). Within that framework, EEG oscillations are the most relevant approach, at least for our research. Although several studies have been conducted to investigate the relationship between EEG oscillations and respiratory signals, it is still unclear how respiration and human alpha activity are related. In a recent review paper, Brændholt et al. (2023) presented a list that included authors, objects (models), frequency, and brain regions. Among all the frequency bands analyzed (LF, delta, theta, alpha, beta1, beta2, and gamma), only three studies refer to human alpha activity: (Perl et al., 2019) EEG alpha and beta, the parietal cortex, and hippocampus; (Kluger and Gross, 2021) MEG delta, theta, alpha, beta, gamma; cortical/subcortical structures, and (Kluger et al., 2021) MEG alpha and parieto-occipital cortex. However, even in studies where the alpha activity was registered, the measured quantity was mostly amplitude or power (Bing-Canar et al., 2016; Herrero et al., 2018; Perl et al., 2019; Kluger et al., 2021). Since their discovery, alpha waves are maybe the most studied cortical oscillations. Current research suggests that they are generated by inhibitory thalamic interconnection, thalamo-cortical feedback loops, and cortico-cortical networks (Da Silva et al., 1980; Nunez et al., 2001). According to new research, cortical structures lead the thalamus, while alpha propagates from higher-order to lower-order areas in both the visual and somatosensory cortex (Halgren et al., 2019). Its amplitude is inversely proportional to the level of cortical activity; strong alpha activity is associated with cortical and behavioral deactivation (e.g., Klimesch, 1999; Klimesch et al., 2007; Rihs et al., 2007). Its involvement in perceptual and memory processes was also studied (Klimesch, 1999; Ergenoglu et al., 2004; Klimesch et al., 2005).

Phases of cortical oscillations play a significant role in brain information processing since they are related to the exact timing of neural activity (Sauseng and Klimesch, 2008). EEG phase synchronization shows communication between distant but functionally related neural populations, where information is exchanged between global and local neuronal networks, as well as the sequential temporal activity of neural processes in response to incoming sensory stimuli (Sauseng and Klimesch, 2008). The relationship between respiratory and EEG alpha phases was studied by Hsu et al. (2020). In their MEG study, the authors explored the effect of slow and normal-paced breathing (0.125 and 0.25 Hz, respectively) on the cosine similarity version of the inter-trial phase coherence (Tallon-Baudry et al., 1996) between respiratory and alpha phases. Based on simultaneous magnetoencephalography and respiratory measurements, they reported that while participants performed slow-paced compared to normal-paced breathing, this condition modulated alpha phase activity, allowing it to appear more organized across a wide range of brain areas. This suggested that slow-paced inspiration was able to organize the cortical alpha phase in a more regularized pattern than the phase associated with normal-paced breathing. Since we found no studies where a respiratory frequency rhythm directly modulates alpha phases (phase differences) in a strictly spontaneous breathing regime, we tried to attempt this research. In our previous papers (Kalauzi et al., 1998; Kalauzi et al., 2009; Kalauzi et al., 2012b; Kalauzi et al., 2018), we proposed, by applying the variable amplitude carrier frequency (CF) mathematical model for human alpha (or rat theta) activity, a procedure to calculate more precisely phase shifts (PS) between any two alpha activity EEG channels than using any Fourier component (FC) phase shift. In this work, for the first time, we applied fast Fourier transform (FFT) on CFPS time courses, FFT (CFPS(t)), in order to examine their oscillatory content.

Body–brain interactions in different states of alertness are recently gaining attention. In addition to the normal resting wake state, these interactions were studied in normal (Abdalbari et al., 2022) and disordered (Faes et al., 2016) sleep. In contrast, a network approach was applied to quantify the cardio–respiratory–EEG relationship in states of heightened alertness: information dynamics (storage and transfer; Zanetti et al., 2019) and multivariate correlation (Pernice et al., 2021). We could not find the analysis of these interactions in the state of drowsiness, a transient state between wakefulness and sleep. Moreover, although several papers address the dependence of human cognition on the respiratory phase (Perl et al., 2019), the topographic dependence of alpha inter-channel phase shift oscillations in different states of alertness was not studied. Therefore, after a detailed validation of the method, topographic distributions of the most prominent amplitudes from this analysis, calculated in wake and drowsy states, are presented and their difference is statistically compared. The state of drowsiness (or hypnagogic state) is divided into nine stages (Hori et al., 1990). Psychophysiological studies have revealed that the main characteristics of this state include decay of alertness, decreased attention, vigilance, and sensory signal detection (Naatanen and Picton, 1987; Ogilvie and Wilkinson, 1984; Ogilvie et al., 1991, respectively). Ogilvie and Wilkinson, (1984) discovered that the transition from wakefulness to sleep is a gradual process. Ogilvie et al. (1991) and Harsh et al. (1994) concluded that the first significant changes in the attention level occurred during Hori’s stages 2 and 3. These data suggest that the first half of the hypnagogic period (Hori’s stages 1–4, where stage 1 has maximal and stage 4 minimal alpha power, without other frequency bands), with its distinct characteristics of decreased attention and dream-like mentation (Nielsen et al., 2005), warrants particular scientific investigation. In our experiments, the state of alertness was monitored by a neurologist who prevented subjects from transitioning to S1 of NREM sleep (N1), i.e., passing beyond Hori’s stage 4. This was confirmed by the fact that the two sets of spectra showed only alpha band amplitude attenuation, without any frequency shift (Bojić et al., 2010).

It should be noted that in this work, we hypothesize that the main rhythm depicted in the respiratory frequency range is caused by respiration. However, since we did not analyze concurrent EEG and respiratory recordings from the same group of individuals, at this stage, we cannot definitely confirm the hypothesis. Therefore, all considerations about body–brain interactions included in this work should be considered conditional, assuming that our hypothesis is correct.

 2 Methods 2.1 Subjects, experimental protocol, and data preparation 2.1.1 EEG recordings

Signals analyzed in this work were originally recorded by Vuckovic et al. (2002) and used in our previous works (Bojić et al., 2010; Kalauzi et al., 2012a; 2012b; 2015; 2018). All recordings were performed in accordance with the medical ethical standards after the subjects signed the informed consent form approved by the local ethical committee. Ten adult healthy human subjects (seven male and three female individuals), age 25–35 (mean ± SD = 28.3 ± 6.5) years, of normal intelligence and without mental disorders were reported after passing a neurological screening. The subjects were lying in a dark room with their eyes closed (standard eyes closed no-task condition; Stam et al., 2002; Stam and van Dijk, 2002). A neurologist monitored their state of alertness and prevented them from falling asleep beyond S1 of NREM sleep. The participants were not previously subjected to sleep deprivation or deviation from their circadian cycles. They were not taking any medication. Based on a visual inspection, artifacts were removed manually. Other details about data collection and preparation can be found in Vuckovic et al. (2002). Two neurologists were independently classifying all signals into wake and drowsy periods. The study included only those sequences that were classified as clearly awake or drowsy by both experts (60 s for each state and each subject).

An EEG machine (MEDELEC 1A97 EEG system, MEDILOG BV, Nieuwkoop, the Netherlands) was used in an electromagnetically shielded room during each 30 min session. EEG electrodes (Ag/AgCl; impedances R _ 5 k_) were positioned at 14 locations (F7, F8, T3, T4, T5, T6, F3, F4, C3, C4, P3, P4, O1, and O2), following the International 10–20 system, with an average reference. The recordings were band-pass filtered between 0.5 and 70 Hz. The EEG recordings were digitized using 12-bit resolution and a sampling rate of 256 Hz per channel (A/D PCI board, Data Translation 2801, Marlboro, MA, USA).

 2.1.2 Respiratory signal recordings

We conducted the experimental protocol on 20 healthy adult human subjects (13 male and 7 female individuals), with mean ± SD = 34.4 ± 7.4. The protocol was approved by the Ethical Committee of the Faculty of Medicine, University of Belgrade (No. 2650/IV-24). The inclusion criteria are as follows: age between 20 and 45 years and absence of any health problems. The exclusion criteria are as follows: patients receiving any medical therapy; a history of pulmonary, cardiovascular, or any other disease at the time of or leading up to the experiments (such as cold, flu, pollen allergy, high temperature, and migraines); and the presence of pathological symptoms during the measurements (high blood pressure, arrhythmias, headache, fatigue, etc.). An additional criterion of exclusion for female participants was the second part of the menstrual cycle (because of its impact on cardiovascular autonomic regulation; Bai et al., 2009; Javorka et al., 2018). All participants were advised to refrain from food and drink for 4 h prior to the start of measurements, to relax and stay attentive, and not to exercise (running, gym, yoga, etc.). Five participants (out of 25) were excluded because of pathological symptoms being discovered during the recordings.

The protocol was executed under controlled laboratory conditions at the Laboratory for Biosignals, Institute for Biophysics, Faculty of Medicine, University of Belgrade. It was conducted between 8 and 12 a.m. in a quiet, refreshing environment at a constant temperature (22°C ± 1°C). All participants were subjected to 10 min of relaxation in a supine position before recording, with no restriction imposed on the air flow. Subjects were advised to adjust the ventilation at the most comfortable rate. They were also instructed not to talk during the recording. Respiration signals were recorded for 20 min. Acquisition was done using the BIOPAC MP100 system (BIOPAC System, Inc., Santa Barbara, CA, USA; AcqKnowledge 3.91 software). The belt with a resistive strain gauge transducer for the continuous recording of breathing was placed slightly above the costal line. All signals were sampled at a frequency rate of 1,000 Hz. We adjusted filters according to Biopac’s instructions for general measurements: gain setting 10, low pass filter with 10 Hz, and no high pass filter (DC-absolute respiratory measurement). Breath-to-breath (BB) intervals were calculated by subtracting the time coordinates of successive inspiration onsets.

 2.2 Data analysis

All programs implementing methods described in this work are original and were developed in MATLAB 2010a (MathWorks Inc., Natick, MA, United States). Statistical tests (Wilcoxon matched-pairs test) were performed using STATISTICA 8.0 (Stat Soft Inc., Tulsa, OK, United States).

2.2.1 A short reminder on the <italic>CFPS</italic> calculation

It is assumed that two oscillations, with equal stable frequencies but variable amplitudes, are being analyzed. The task of determining their phase shift, regardless of whether it be at a single point in time or during the whole recorded time course, is hampered by the fact that each of the two oscillations is spread across a frequency region. The question as to which FC phase shift to rely on is raised. As described in detail in our previous works (Kalauzi et al., 1998; Kalauzi et al., 2009; Kalauzi et al., 2012b; Kalauzi et al., 2018), there are two reasons why a single oscillatory activity, with a stable frequency but variable amplitude, occupies a range of frequencies (rather than only one) when subjected to a concrete FFT analysis. The first reason is that its frequency is positioned between two adjacent FCs’ frequencies (which is almost always the case) and which we conveniently named “inter-componentality.” The second reason originates from the variability of amplitudes, where this variability raises to two sidelobes in the amplitude spectrum (in radio communications termed double sideband amplitude modulation, DSBAM).

For the first cause, let us observe a cosine wave with constant frequency f c = ω c /(2π) (usually termed carrier frequency, CF), constant amplitude A c , and initial phase φ c : u t = A c ⁡ cos ω c t + φ c .

If such a signal of duration T is subjected to FFT analysis, then CF can be expressed as ω c = k ω 0 + Δ ω , where ω 0 = 2π/T and kω 0 denotes the nearest lower Fourier component frequency. The amplitude spectrum of the signal u(t) can then be expressed as A m n ω 0 = x n 2 + y n 2 , where real and imaginary parts of the spectrum are x n = 4 A T k ω 0 + Δ ω k ω 0 + Δ ω 2 + n ω 0 2 cos φ c + π Δ ω ω 0 sin π Δ ω ω 0 , y n = − 4 A T n ω 0 k ω 0 + Δ ω 2 + n ω 0 2 sin φ c + π Δ ω ω 0 sin π Δ ω ω 0 .

It is clear from Equations 1–3 that this single-frequency oscillation occupies the whole available frequency range, with the largest amplitudes occurring next to ω c and decreasing as one moves away from it.

Subsequently, let us observe two such oscillations with equal carrier frequencies: u c 1 t = A c 1 ⁡ cos ω c t + φ c 1 ;   u c 2 t = A c 2 ⁡ cos ω c t + φ c 2 .

Equations 1–3 yield the following phase shift spectrum: Δ φ n ω 0 = tan − 1 n ω 0 S 1 ω c C 2 − ω c C 1 n ω 0 S 2 ω c C 1 ω c C 2 + n ω 0 S 1 n ω 0 S 2 , where S i = sin φ c i + π Δ ω ω 0 ;   C i = cos φ c i + π Δ ω ω 0 ; i = 1 , 2 .

Based on Equation 4, the expansion of Δφ(nω 0 ) into Taylor series around ω c shows that there is an approximately linear dependence of Δφ(nω 0 ) on nω 0 in the vicinity of ω c , with FC phase shift values being around the carrier frequency phase shift Δφ c = φ c2 – φ c1 (Kalauzi et al., 2012b, Appendix B). “Being around” means that Δφ(nω 0 ) – Δφ c changes sign when nω 0 – ω c does. Therefore, by angularly averaging FC phase shifts around ω c , we can obtain a more accurate estimate of CFPS, or Δφ c , than by relying on any of the single FC phase shifts.

Regarding the second cause, i.e., the influence of amplitude variability on FC phase shifts, let us again observe two CF signals, this time both with variable amplitudes, where amplitude variability is achieved by applying amplitude modulation on both signals: u c i t = 1 + u m i t U c i ⁡ cos ω c t + φ c i , i = 1 , 2 , where u mi stands for the two modulating functions. After expanding both modulating functions into Fourier series and keeping the first n FCs, Equation 5 can be approximated with u c i t ≈ 1 + ∑ j = 1 n U m j , i ⁡ cos ω m j t + φ m j , i U c i ⁡ cos ω c t + φ c i , i = 1 , 2 , where U mj,i , ω mj , and φ mj,i stand for the amplitude, frequency, and initial phase of the jth FC of the ith modulating function, respectively (please note that the modulating FC frequency ω mj does not depend on i).

After multiplication, if we transform cosine products into sums, we obtain u c i t ≈ U c i ⁡ cos ω c t + φ c i + 1 2 ∑ j = 1 n U m j , i U c i [ cos ω c + ω m j t + φ c i + φ m j , i + cos ω c − ω m j t + φ c i − φ m j , i ] .

Notably, Equation 6 contains three terms: one component at CF, ω c , and two sidebands at ω c + ω mj and ω c –ω mj . Components at the two sidebands have initial phases φ ci + φ mj,i and φ ci –φ mj,i , with an opposite sign of the second term. Therefore, the CF initial phase, φ ci , can again be obtained by angularly averaging phases of FCs around the carrier frequency, ω c . The same is valid for two signals, where by angular averaging of FC phase shifts around CF, we can calculate CFPS more accurately than by relying on any of FC phase shifts. Finally, according to both Equations 4, 6, in order to calculate CFPS of two signals as accurately as possible, both causes of the signals occupying a frequency range can be overcome by the same procedure–angular averaging of FC phase shifts around the CF. Accuracy of this calculation is further increased if we perform a weighted averaging with FC amplitudes as weights.

 2.2.2 Angle range extension method

In our previous works, we treated CFPS as a random variable, working mostly with distributions and histograms and deriving basic statistical properties from them (Kalauzi et al., 2009; Kalauzi et al., 2012b). Exploring time dependence of this quantity, CFPS(t), and its FFT spectral properties was not studied until now. In order to do a successful FFT analysis, we must first introduce some new procedures for preprocessing CFPS(t). One important step in that direction is described in the following section.

In this work, we performed two FFT procedures: the first procedure (“primary”) on initial EEG recordings. Since the duration of each recording was 60 s, with a primary sampling frequency of 256/s, 15360 samples were generated. They were subjected to a moving window FFT analysis, using 2-s long overlapping moving windows that advanced at a 0.25 s step. With these parameters, for each subject and each EEG channel, 233 spectra could be obtained, producing CFPS(t) waveforms with the same number of points. Secondary FFT analysis was done on these waveforms, but because they were relatively short, each of the 910 CFPS(t) signals was analyzed as a whole, without any moving windows, but after previous detrending. Although we expected relatively slow CFPS(t) oscillations to appear, we chose overlapping windows in order to generate a sufficiently high secondary sampling frequency (4 Hz), allowing us to detect CFPS(t) oscillations of up to 2 Hz. Analysis parameters of both levels (primary and secondary) are summarized in Table 1.

Summary of parameter values for the two levels of analysis.

Primary Secondary
Input signal EEG CFPS(t)
Output signal CFPS(t) FFT (CFPS(t))
Sampling frequency 256/s 4/s
Window length (FFT epoch) 2 s (60–2)/0.25 + 1 = 233 spectra233*0.25 = 58.25 s
Window step 0.25 s No moving window, whole signal FFT

When calculating formulas involving angles, it is often the case that the results are limited to the range of [-180°, 180°] or, alternatively, [0°, 360°]. An example of this fact is the often applied MATLAB atan2 command. In some physiological experiments dealing with oscillations, relatively small angle deviations around the zero baseline appear, requiring no correction of angle range. However, in many other cases, deviations of phases or phase-shifts of EEG signals are substantial. If time courses of angles are being recorded, a considerably artificial limitation to the abovementioned ranges causes frequent and big, abrupt jumps of angle values, as shown in angle time plots (Figure 1A). In order to achieve a more natural time dependence of angles, a correction procedure, named angle range extension (ARE), is proposed for the first time. This procedure is essentially aimed to detect these “abrupt jumps”, i.e., differences between two adjacent CFPS(t) points, which are above a certain critical value (we named it “critical phase shift difference”, C psd ). If CFPS(i) is the CF phase shift calculated at the ith Fourier window position, then if CFPS(i+1) < CFPS(i) and abs (CFPS(i+1) - CFPS(i)) > C psd , then CFPSc(i+1) = CFPS(i+1) +360°, where CFPSc(i+1) stands for the corrected value. Alternatively, if CFPS(i+1) > CFPS(i) and abs (CFPS(i+1) - CFPS(i)) > C psd , then CFPSc(i+1) = CFPS(i+1) −360°. Therefore, addition is performed in case of a sudden decrease and subtraction in case of an increase. However, determining the optimal value of C psd remains a challenge. It can be assumed that the corrected CFPS(i) waveform should have a smaller line length, ∑ i abs C F P S i + 1 − C F P S i , than the uncorrected line since long jumps, supposedly, were eliminated during the procedure. Therefore, for each waveform, we applied a series of ARE corrections, where C psd was varied. For each C psd , the corrected line length was calculated, and C psd , where the line had a minimal length, was taken as the optimal value (Figure 2).

Representative example of the ARE method, in which the procedure, applied on the CFPS(t) waveform, had a substantial effect. (A) Uncorrected, derived from subject one between channels T3 and O2, where a contra-phase was present (CFPS values grouped approximately ±180°). (B) Corrected after applying the ARE procedure, where most of the “abrupt jumps” were eliminated, making it more suitable for further FFT analysis. The contra-phase is present usually when fronto-occipital or temporo-occipital pairs are being analyzed.

Example of dependence of the corrected CFPS line length on the critical phase shift difference, C psd . Optimal C psd , above which the ARE procedure is to be applied, corresponds to the minimal CFPS line length and is marked with a red dot. This example refers to the CFPS waveform presented in Figure 1A to produce the corrected waveform shown in Figure 1B.

 2.2.3 FFT analysis of carrier frequency phase shifts, <italic>CFPS</italic>(<italic>t</italic>)

In our previous work (Kalauzi et al., 2012b), we studied statistical properties of CFPS of 10 test individuals in the wake and drowsy states, treating it as a random variable. In this study, we aimed at investigating probability distributions of phase shifts or phase potentials, from which properties such as the average value, most probable value, and standard deviation could be derived. At that stage, we did not analyze oscillatory properties of CFPS(t) waveforms, which could yield information about possible links of EEG with other physiological rhythms, such as ECG and respiration (so-called “mind–body” relationship).

Currently, we tried to combine two methods:

a) By treating alpha EEG activity as a carrier-stable inter-componential frequency and amplitude variable signal, we perform CFPS calculation, integrating Fourier component phase shifts from the whole alpha frequency range (8–12 Hz), according to the procedure described in our previous work (Kalauzi et al., 2012b).

b) On these “raw” CFPS waveforms, we apply the ARE method in order to obtain the waveforms suitable for further FFT analysis (parameters shown in Table 1).

2.2.4 Method validation

In this work, we tried to validate our method by subjecting an artificial, mathematically generated signal to the same procedure as the recorded CFPS(t) time dependencies. To achieve this, we first simulated the profiles of relative FFT amplitudes of typical CFPS(t) spectra. Two of them are presented in Figure 3.

Two representative examples of relative amplitudes of CFPS(t) FFT spectra. Profiles are similar, exhibiting different individual oscillations superimposed on an fGn + fBm type noisy background. (A) Channel pair F8, T4 recorded in subject 7. (B) Channel pair F8, O2 is recorded in subject 9.

As can be noticed, spectra are complex, consisting of different individual oscillations superimposed on a noisy background. However, these individual peaks did not appear in all 910 spectra. Specifically, they were more visible in spectra of some of the channel pairs than those in others, while altogether missing in some. In addition, inter-individual variability prevented them from appearing in all subjects for a given pair of EEG channels. Despite this complex behavior, in order to validate the method, it is necessary to perform the following five steps:

a) Generate an appropriate simulated background (of the corresponding fractal nature), as well as one individual artificial oscillation, preferably in the frequency range 0 < f < 0.5 Hz (according to the profiles shown in Figure 3), and add them mathematically.

b) Generate two synthetic 10-Hz carrier frequency “alpha” signals, with arbitrary constant (or random) amplitudes (Bojić et al., 2010).

c) Perform phase modulation of one of the signals in b) with the signal described in a).

d) Perform phase demodulation, i.e., calculate CFPS(t) (according to the procedure described in our previous works, Kalauzi et al., 1998; Kalauzi et al., 2012b) between the modulated and unmodulated synthetic “alpha” signals, perform the ARE procedure, and calculate the FFT spectrum of the resulting CFPS(t) time series.

e) Compare and statistically test the original signal in a) against the demodulated CFPS(t), as well as their FFT spectra.

2.2.4.1 Simulating the <italic>CFPS</italic>(<italic>t</italic>) background activity

In order to identify to which type of signals, such as two CFPS(t) examples shown in Figure 1, belong to, we reverted to the well-known but simple criteria explained in some previous papers (Eke et al., 2000) and employed in some of our previous works (Bojić et al., 2010), where the value of the exponent β determines its fractal nature. This exponent is given by the linear regression slope of the signals’ power vs. the frequency log–log plot: A f 2 ∝ c f − β , where c stands for the proportionality factor. According to this classification, the critical value of β is 1 since for −1 < β < 1, the signal is identified as the fractal Gaussian noise (fGn), which is a stationary signal with constant variance. On the other hand, if 1 < β < 3, the time series is regarded as the non-stationary fractal Brown motion (fBm), in which signal variance increases with its length: v a r X t ∝ t 2 H , where X denotes the current signal value and H is the Hurst exponent (Mandelbrot and van Ness, 1968; Eke et al., 2000).

As we could not find any specific data about the fractal nature of CFPS(t) in the literature, we performed linear regressions of all 1820 log–log power vs. frequency CFPS(t) signals (91 channel pairs × 10 subjects × 2 states) and histogrammized the corresponding beta values. The results are presented in Figure 4.

(A) Histogram of all 910 β-values (91 channel pairs × 10 subjects) calculated according to Equation 7. (B) Histogram of 91 β-values calculated for each channel pair by averaging the values across the group of 10 subjects, according to Equation 8. Vertical red lines mark the mean β.

Let β s , s t c 1 , c 2 denote this quantity between channels c1 and c2 of subject s in state st = {wake, drowsy}, which was obtained from CFPS(t) using Equation 7. As shown in Figure 4A, all 910 such values were histogrammized, showing a unimodal distribution around mean(β) = 1.6596 ± 0.2741. Interestingly, the majority (98.57%) of the values were 1 < β < 2.5 and could be characterized as fBm, with only 13 (1.43%) belonging to the fGn type (β < 1). However, we average individual beta values for each channel pair across the tested group and denote them as β s t c 1 , c 2 = 1 10 ∑ s = 1 10 β s , s t c 1 , c 2 . The corresponding histogram (Figure 4B) points to a bimodal distribution: one mode localized at approximately β 1 = 1.59 and the other at approximately β 2 = 1.72. In this case, the overall mean is (β) = 1.6741 ± 0.1242. However, all these β-values, although not originating from individual signals but were a result of averaging, had an “fBm nature” (1 < β < 2).

There are different ways to generate a synthetic signal with desired fractal properties, such as using Weierstrass functions (Esteller et al., 2001; Kalauzi et al., 2005), in which the signal fractal dimension (D) is theoretically given and directly linked to one of the two parameters (γ, H): W H γ t = ∑ t γ − i H ⁡ cos   2 π γ i t , where 0 < H < 1, γ > 1, and H = 2 - D.

However, the range of β, found in the histogram abscissa of Figure 4A, suggests that both fGn and fBm contribute to the FFT background of CFPS(t) signals. This circumstance may point to an alternative way of generating a signal with the desired β-value—a composite fGn/fBm. This fractal signal can be obtained by summing the two componential signals with various relative amplitude contributions. For this purpose, we used the approach described in Kalauzi et al. (2012a), where two signals with inverse amplitudes were summed while keeping the amplitude of the sum constant: A f G n f B m = 1 − k a A f G n + k a A f B m , where A fGnfBm stands for the amplitude of the composite fGn + fBm signal, A fGn and A fBm denote amplitudes of the two componential signals, and k a is the relative amplitude attenuation factor (0 ≤ k a ≤ 1). A custom MATLAB program was made to calculate how β of the composite fGn + fBm signal depends on k a . The result is presented in Figure 5. In order to cover the whole beta range (0 < β < 2), fGn was represented by a synthetic white noise Wn(i), i = 1,.,233, acting as a random variable with the uniform probability distribution U (−0.5, 0.5). Then, we derived the Brown motion signal, Bm(j), j = 1,.,233, representing the fBm component, by numerically summing white noise samples: B m j = ∑ i = 1 j W n i , j = 1,.,233.

(A) Dependence of exponent β (Equation 7) on the relative amplitude attenuation factor k a (Equation 9) for 20 repeatedly generated white noise and 20 Brownian motion signals. Thin lines represent 400 of these repetitions, whereas the thick black line represents their average. (B) Average (solid) ± standard deviation (dashed) of the same set of lines as shown in panel (A). Mean β (1.6596), from the histogram shown in Figure 4A, is plotted as a thin horizontal line in order to determine the corresponding k a value (≈0.92) to be used in the CFPS(t) background simulation.

Subsequently, we aimed to establish a β = f (k a ) dependence, from which, using it as a sort of calibration line, the value of k a , corresponding to the mean β-value, obtained from the 910 CFPS(t) signals (mean(β) = 1.6596; Figure 4A), could be determined. Since both Wn and Bm signals, generated with a relatively small number of samples in order to match the recorded signals, might show some statistical variability, we repeated the construction of the β = f(k a ) calibration line using 20 Wn and 20 Bm signals, obtaining a field of 400 combination lines (Figure 5A). In addition, since these lines also depend on the initial amplitudes of both Wn and Bm components, amplitude normalization had to be carried out: y N i = y i − m e a n y s t d y , where y(i) stands for the ith sample value of either of the fractal signals (Wn or Bm).

The importance of precisely determining the coefficient k a for generating a successful background simulation is illustrated in Figure 6. Figure 6A shows an example of the recorded CFPS(t) (subject 1; EEG channels F8, T4), while Figures 6B, C contain simulated background activity for k a = 0.92 and k a = 0.5, respectively. The similarity of signals shown in Figures 6A, B is obvious, which is contrary to Figures 6A, C.

Comparison between recorded and simulated CFPS(t) background activities in the time domain. (A) Recorded from subject 1: channels F8, T4. (B) Simulation containing Wn and Bm, with k a = 0.92. (C) Simulation with k a = 0.5. (D) Relative amplitude FFT spectra of 400 repeated simulations generated using k a = 0.92. Thin colored lines represent individual simulations, whereas the thick black line represents their average.

Finally, to the background presented in Figure 6B, a mathematically synthesized sinusoid of 0.22 Hz frequency was added in order to simulate the situation where both the background and the most prominent recorded oscillation was present (such as those shown in Figures 3, 7A). Relative FFT spectra of the simulated, as well as an example of one of the recorded signals, are shown in Figure 7. It is important to note that since the recorded CFPS(t) signals differed among themselves, it was not necessary to simulate any of them exactly; for method validation, it was sufficient that the simulation reproduced signal’s main properties.

(A) Representative example of a relative amplitude spectrum of recorded CFPS(t) activity, with pronounced individual oscillations superimposed on the background, obtained from subject 9, between channels O1 and O2. (B) Simulation of the activity shown in A, consisting of a fGn + fBm (represented by Wn + Bm) fractal background, with k a = 0.92 and summed with a 0.22-Hz artificial sinewave.

Once generated, the simulated signal was used to follow the procedures described in steps a)–e) in Section 2.4 and also used as the input to the CFPS(t) calculation described in detail in our previous papers (Kalauzi et al., 1998; Kalauzi et al., 2009; Kalauzi et al., 2012b) and in this study in Section 2.2. For the method to be validated, the input signal should resemble, as much as possible, the output signal obtained after completing the abovementioned steps. Their comparison is shown in Figure 8. The coefficient of linear correlation between the input and output signals was r = 0.9421, and it was highly significant (p = 1.46 × 10–111). The resulting (output) signal shows (Figures 8A, B) that some higher frequency oscillations were attenuated due to the low-pass filtering effect of the 2-s moving Fourier window (Kalauzi et al., 2023a). However, this attenuation has little effect on our main finding regarding the presence of CFPS(t) oscillations because they are positioned at ≈ 0.22 Hz, while the first filter zero is at 0.5 Hz. Moreover, this attenuation is equal for all 91 channel pairs so that CFPS(t) waveforms and their FFT amplitudes could be directly compared.

(A) Simulation of a typical CFPS(t) waveform, containing two fractal background components, Wn and Bm, and a 0.22-Hz sinewave. The simulation served as the input to the procedure of calculating CFPS(t), which is part of the steps a)–e) in Section 2.4. (B) Output from the CFPS procedure, having a high correlation with the input (r = 0.9421). Inlet: input–output waveforms superimposed.

If we denote with A s , s t c 1 , c 2 f the FFT amplitude spectra of CFPS(t) waveforms between channels c1 and c2, recorded in subject s in state st (analogously to the notation of Equation 8), and average them across the tested group and all channel pairs in order to get one resulting spectrum, we obtain A r s t f ¯ ¯ = 1 910 ∑ c 1 , c 2 = 1 91 ∑ s = 1 10 A r s , s t c 1 , c 2 f , where A r stands for the relative amplitude: A r f = A f ∑ f A f .

Finally, spectra averaging across the group for each channel pair result in the following expression: A r s t c 1 , c 2 f ¯ = 1 10 ∑ s = 1 10 A r s , s t c 1 , c 2 f .

A custom MATLAB program was made to select the given number of channel pairs with spectra having the maximal sum of relative amplitudes in the chosen frequency region.

 3 Results 3.1 Results of FFT analysis applied on ARE-corrected alpha <italic>CFPS</italic>(<italic>t</italic>) waveforms of a test group in the wake state

The spectrum defined in Equation 10 is presented in Figure 9A. The profile of this spectrum shows some particular oscillations superimposed on the formerly analyzed background activity, consisting of mixed fGn and fBm fractal signals, with the corresponding β exponents’ histogram presented in Figure 4A. Amplitude variability visible at f > 0.5 Hz belongs to sidelobes, which result from the low-pass filtering effect of the 2-s FFT moving window. The most visible individual peak is the one positioned at 0.2232 Hz (13th FC). However, individual spectra (from the pool of 910 available) differ among each other in amplitudes and exact positions of these peaks, which appear in the 0 < f < 0.5 Hz region. In order to explore which channel pairs exhibit most prominent relative FFT amplitudes on f = 0.2232 Hz and hopefully increase the signal-to-noise ratio (S/N), the first step was to perform spectra averaging across the group for each channel pair according to Equation 11. If we choose to average and present five channel pairing with the maximal sum of relative amplitudes in the region of 12th–15th FCs (0.2060–0.2575 Hz), according to Equation 11, the resulting spectrum is presented in Figure 9B, where the improved S/N ratio is apparent compared to that shown in Figure 9A. The list of these five channel pairs, as well as the corresponding values of A r w c 1 , c 2 0.2060   –   0.2575 ¯ , is listed in Table 2, while the topographic distribution of these channel pairs is plotted in Figure 10A, where mostly longitudinal (fronto-occipital) channels are connected. Analogous calculation was carried out for one FC (13th) by applying A r w c 1 , c 2 0.2232 ¯ . The result is also given in Table 2 and shown in Figure 10B. As can be seen, pairs of connected EEG channels were similar in both cases.

(A) Average of all 910 (10 subjects x 91 channel pairs) FFT CFPS(t) spectra according to Equation 10; no selection applied. A peak at the 13th FC (0.2232 Hz) is small but visible. (B) Average of five spectra, previously averaged across 10 subjects (Equation 11), selected to have maximal sum of relative amplitudes between 12th and 15th FCs (0.2060–0.2575 Hz). The same peak appeared, with an increased S/N ratio compared to the ratio in panel (A). (C) Average of five spectra, with no previous group averaging (i.e., directly from the pool of 910 spectra), again characterized by the maximal sum of relative amplitudes between 12th and 15th FCs. Main peaks were in the range of 0.2060–0.2403 Hz, as well as one appearing at sixth FC (0.1030 Hz). Inlet: an example of recorded CFPS(t), subject 9, channels O1 and O2, where respiratory frequency oscillation could be noticed by the naked eye. (D) Histogram of 5450 BB (breath-to-breath) frequencies recorded from 20 healthy subjects in the supine position. Vertical red lines, with mean ± st. dev., also drawn in panel (C). Thin colored lines in panels (B, C) originate from five individual spectra, whereas the thick black line represents their average.

Upper part: Five channel pairs with maximal sum of 12th–15th FC relative FFT amplitudes and their corresponding A r s t c 1 , c 2 f ¯ values (group mean ± st. dev.). Lower part: The same calculation for only one (13th) FC.

Wake Drowsy
Channel A r w c 1 , c 2 0.2060   –   0.2575 ¯ Channel A r d c 1 , c 2 0.2060   –   0.2575 ¯
F8, 02 0.1214 ± 0.0396 F8, F4 0.0913 ± 0.0249
F4, 02 0.1187 ± 0.0486 C3, 02 0.0911 ± 0.0223
01, 02 0.1185 ± 0.0503 01, 02 0.0898 ± 0.0372
F3, 01 0.1154 ± 0.0441 C4, P4 0.0872 ± 0.0300
F3, 02 0.1144 ± 0.0354 T5, 01 0.0850 ± 0.0156
Channel A r w c 1 , c 2 0.2232 ¯ Channel A r d c 1 , c 2 0.2232 ¯
F8, 02 0.0427 ± 0.0123 F8, F4 0.0285 ± 0.0132
F8, T6 0.0373 ± 0.0194 C3, 02 0.0280 ± 0.0091
F4, 02 0.0367 ± 0.0170 01, 02 0.0279 ± 0.0145
F3, 02 0.0348 ± 0.0157 C4, P4 0.0257 ± 0.0106
P4, 02 0.0332 ± 0.0140 F8, F3 0.0252 ± 0.0131

Topographic view of the selection of five channel pairs with a maximal group average of relative FFT amplitudes in the wake state, according to Equation 11. (A) Channel pairs with the maximal sum of four relative amplitudes (12th–15th FCs), corresponding to the frequency range of 0.2060–0.2575 Hz. (B) Analogous topography shown in panel A drawn for one (13th FC) at a frequency of 0.2232 Hz. (C) Relative FFT amplitudes for all 91 channel pairs (Equation 11) plotted in the descending order. Ordinates of the red line (representing one FC) were taken from Table 2 and multiplied by four to make them comparable to the blue line (corresponding to summed amplitudes of 4 FCs). The first five channel pairs of the red line are marked with black circles because they correspond to the pairs plotted in panel (B). Their abrupt decrease explains why we chose five pairs for the topographic presentation. For the sake of clarity, error bars were omitted.

Finally, if we want to eliminate the influence of inter-individual variability on oscillation detection and appearance, and perhaps even further improve the S/N ratio, no primary averaging should be performed, and five spectra with the maximum sum of relative amplitudes in the 0.2060–0.2575 Hz region should be selected directly from the pool of all 910 available spectra. The resulting average of these five finally selected spectra is presented in Figure 9C. In this case, main peaks were in the range of 0.2060–0.2403 Hz; however, a peak at the sixth FC (0.1030 Hz) also appeared. A further increase in the S/N ratio was herewith achieved, with even less fractal-type background activity present. In order to identify, not conclusively but to a reasonable degree of probability, the possible physiological origin of the oscillation positioned at ≈ 0.21–0.24 Hz, we used breath-to-breath (BB) interval data from our previous experiments (Matić et al., 2020; Matić et al., 2022; Kalauzi et al., 2023b), where respiratory signals from 20 healthy subjects in the supine position were recorded. A histogram of their 5,450 pooled instantaneous breathing frequencies (B f ), which were obtained as the inverse of the corresponding BB intervals, is presented in Figure 9D. These pooled data had a mean (B f ) = 0.2552 ± 0.0715 Hz and are indicated as vertical red lines in Figures 9C, D. From this comparison, it can be seen that the position of the ≈0.21–0.24 Hz FFT CFPS(t) peak lies within the expected frequency boundaries of relaxed healthy human breathing, although slightly shifted toward lower frequencies. If the breathing origin hypothesis is adopted, the small shift could be explained by the fact that the 10 EEG subjects were more relaxed than the 20 respiratory signal subjects (Matić et al., 2020; Matić et al., 2022) as the former were in the wake state but not far from the drowsy state (Vuckovic et al., 2002; Kalauzi et al., 2012b). Final confirmation of the respiratory connection of detected CFPS(t) oscillations is expected to be obtained by simultaneously recording both respiratory and EEG signals from a sufficient number of subjects, preferably in different cardio-respiratory conditions, where frequency matching and phase locking between respiratory and CFPS(t) oscillations can be accurately established.

In addition to topographic distributions on Figures 10A‐C, show group averages of relative FFT amplitudes for all 91 channel pairs, arranged in descending order using Equation 11. Ordinates of the red line (referring to 1 FC) were derived from Table 2 and multiplied by 4 to make them comparable to the blue line, which depicts the group average of four summed FC (12th–15th) amplitudes. Notably, first five channel pairs of the red line are marked with black circles because they correspond to the pairs plotted in Figure 10B. Five pairs were chosen because they belong to the part of the line where an abrupt decrease occurs.

 3.2 Results of FFT analysis applied on ARE-corrected alpha <italic>CFPS</italic>(<italic>t</italic>) waveforms of a test group in the drowsy state

The same procedure was applied to the recorded signals from the same 10 subjects in the drowsy and wake states, allowing the results to be directly compared.

Drowsy state results are displayed in Figures 11, 12. Figures 11A–C essentially show the same information as that presented in Figures 9A–C. Figure 11D shows a direct comparison of averaged spectra from those shown in Figures 9C, 11C. From all panels, as shown in Figure 11, one can observe a reduction in individual peaks in the low frequency region (f < 0.3 Hz) compared to the wake peaks shown in Figure 9.

(A–C) Spectra corresponding to those shown in Figures 9A–C, but here, they refer to the drowsy state. (A) CFPS oscillation peaks in the low-frequency region f < 0.5 Hz are reduced or even diminished compared to the wake state, as shown in Figure 9A. (B) After averaging across the 10 subjects, in five selected spectra, a reduced 13th FC (0.2232 Hz) peak appeared in the drowsy state as well. (C) The sixth FC peak at 0.1030 Hz, which was present in the wake state shown in Figure 9C, disappeared in the drowsy state. (D) Averaged spectra from Figures 9C, 11C are superimposed so that they can be directly compared.

(A, B) Same topographic distribution of five channel pairs is shown in Figures 10A, B, obtained for the drowsy state. As shown in the wake state, the two distributions are mutually similar but differ significantly from the distributions in the wake state. (C) Relative FFT amplitudes for all 91 channel pairs (Equation 11) plotted in the descending order, superimposed for the wake (w) and drowsy (d) states and for two frequency regions (13th and 12th – 15th FCs). As shown in Figure 10C, ordinates of the brown line (one FC) were multiplied by four to be comparable to the black line (summed amplitudes of 4 FCs). For all channel pairs, amplitudes in the drowsy state are smaller than those in the wake state, pointing to the fact that wakefulness is necessary for the respiratory frequency alpha CFPS to be present in full capacity across the whole scalp.

The topographic distribution of channel pairs, characterized by five greatest group averages of summed four relative FFT amplitudes in the range of 0.2060–0.2575 Hz (12th to 15th FCs; Figure 12A), or one FC positioned at 0.2232 Hz (13th FC; Figure 12B), was mutually similar but differed significantly from the two distributions obtained in the wake state (Figures 10A, B). Specifically, pairs connecting fronto-occipital regions, recorded in the wake state, did not emerge in the drowsy state.

 3.3 Statistical testing of relative amplitude differences between the wake and drowsy states

The channel pairs shown in Figure 12C were not arranged and presented in the same order since each line reflected its own descending order, rendering them unsuitable for statistical testing. Therefore, we rearranged the order of channel pairs of the drowsy state (corresponding to the order on the brown line shown in Figure 12C) to match the descending wake order (corresponding to the red line on the same panel). Using the Wilcoxon matched-pairs test, we tested differences between wake and drowsy group ensembles of FFT (CFPS(t)) relative amplitudes at 0.2232 Hz (13th FC) for each of the 91 channel pairs. Therefore, each of the two input ensembles for the 91 tests consisted of 10 values, one for each subject, defined in Section 3.1 as A r s , w c 1 , c 2 f and A r s , d c 1 , c 2 f , with s = 1,.,10. Here, c1 and c2 stand for the two EEG channels, while w and d subscripts denote wake and drowsy states, respectively. Of the 91 tested pairs, 20 appeared to have significantly greater (p < 0.05) amplitudes in the wake state. All group ensemble values, averages of relative amplitudes in wake and drowsy states, and the corresponding p-values can be obtained upon request from the authors. In this study, we shall present only five pairs with the most significant differences (p < 0.01): F402 and T4F4, p = 0.0051; T3T6 and F7F8, p = 0.0069; and F8T6, p = 0.0093. The topographic distribution of all 20 pairs is presented in Figure 13. As can be seen, diagonal directions were preferred, with more lines connecting left anterior to right posterior scalp regions than those in the opposite direction.

Topographic distribution of channel pairs which had significantly greater (p < 0.05) relative FFT amplitudes of CFPS(t) oscillations at 0.2232 Hz in the wake state than those in the drowsy state.

 4 Discussion

In this work, the carrier frequency mathematical model of human EEG alpha oscillations was applied to 10 young healthy participants in wake and drowsy states to calculate the time course of CF phase shifts CFPS(t) for all 91 channel pairs. After the validation of the method, their oscillatory content was studied through FFT (CFPS(t)). Although not appearing equally in all channel pairs and every subject, a clear peak in the respiratory frequency region, 0.21–0.26 Hz, was observed (max at ≈ 0.22 Hz). When five channel pairs with the most prominent group averaged amplitudes at 0.22 Hz were plotted in both states, topographic distributions changed significantly—from longitudinal, connecting frontal and posterior channels in the wake state to topographically split two separate regions—frontal and posterior in the drowsy state. Moreover, in the drowsy state, 0.22-Hz amplitudes decreased for all pairs, while statistically significant reduction was obtained for 20/91 (22%) pairs.

4.1 Alpha waves and neuronal excitability

Both the alpha amplitude and instantaneous phase influence neuronal excitability in the cortex. Sauseng et al. (2009) found that a motor-evoked potential was elicited more easily when alpha power immediately preceding the transcranial magnetic stimulation (TMS) pulse was low and vice versa. Similarly, Romei et al. (2008) reported that low pre-stimulus alpha-band power resulted in TMS reliably inducing positive responses, whereas high pre-stimulus values caused the same TMS stimuli to fail to evoke a visual percept during spontaneous fluctuations in posterior alpha activity. However, there are some contradicting results concerning the pre-stimulus amplitude (but not the phase). Berger et al. (2014) applied single-pulse transcranial magnetic stimulation (TMS) over the left sensorimotor cortex while simultaneously recording EEG. Results indicate that the instantaneous phase, but not amplitude, of oscillations in various frequency bands at the stimulation site during the time of TMS pulse correlates with different levels of excitability. The fact that neuronal excitability was dependent on the instantaneous alpha phase had been recognized for quite some time. Callaway and Yeager (1960) found that visual reaction time depended on different phases of the alpha cycle. More recently, Mathewson et al. (2009) showed that both subsequent visual detection and stimulus-elicited cortical activation levels in a metacontrast masking paradigm can be reliably predicted by the phase of alpha rhythm measured over posterior brain regions. When a target was presented during the trough (minimum) of an alpha wave, cortical activation was suppressed 100 m after the presentation, and participants were less likely to detect the target. They concluded that during one alpha cycle, the human brain goes through rapid oscillatory shifts in excitability, directly influencing whether an environmental stimulus will reach conscious awareness or not.

 4.2 Respiration and neuronal excitability

On the other hand, respiration exerts its own impact on neuronal excitability. Neuronal activity phase-locked to respiration has been evidenced in many published results. Respiratory modulation of the power of gamma frequency oscillations was found in the whisker barrel cortex of wake mice (Ito et al., 2014; Heck et al., 2017). Jung F. et al. (2023) stated that “rhythmic body processes strongly influence activity patterns throughout the brain,” further mentioning respiration as a major factor. This influence is extended to both local field potential and single-cell levels. However, according to the authors, it is not known how respiratory-driven rhythms interact or compete with internal brain oscillations, especially if the frequency ranges overlap, such as respiratory and theta rhythms in mice. They report state-dependence of parietal cortex unit spike modulation by both rhythms: during REM sleep, theta dominated unit discharge, while respiration exerted greater influence during active waking. Heck et al. (2017) argued that “respiration, via multiple sensory pathways, contributes a rhythmic component to the ongoing cortical activity,” suggesting that “this rhythmic activity modulates the temporal organization of cortical neurodynamics, thereby linking higher cortical functions to the process of breathing.” Let us note that measurements of respiratory influence on the power of field potential oscillations or unit discharges, such as those in these three studies, differ in one important aspect from the approach involving phases, obtained by, e.g., Hsu et al. (2020) or ours, presented in this work. In other words, in the former case, each particular result originates from one static electrode, resulting in “static” results on the receiving end. In the phase (shift) approach, the results are “dynamic,” describing the interaction of phases of two oscillations: respiratory and field potential.

The combined results of these two lines of research indicate that in this work, we are dealing with two very complex interacting oscillatory systems, each one altering the excitability of extensive neuronal networks and affecting simultaneously all cortical areas. At this stage, when the phase–frequency interaction of the two excitability altering oscillations has only recently been detected, it is not easy to give a valid interpretation or give a meaningful proposal for either anatomical or physiological mechanisms that underlie the phenomenon.

However, some aspects of the issue can still be discussed even in this preliminary phase. Logically, minimal anatomical and consequently functionally active set of structures should include the PreBötzinger complex (preBötC) locus coeruleus and thalamo-cortical system. Influences between respiratory centers and the cerebral cortex are known to be bilateral. PreBötC projects respiratory afferents to cortical brain areas, e.g., via the locus coeruleus, from olfactory nuclei and through the vagus nerve (Del Negro et al., 2018; Kluger and Gross, 2020). In addition, the inspiratory rhythm originating from the preBötC activity coordinates associated conditional oscillators for other phases of the breathing cycle, for sighing and orofacial behaviors. Furthermore, the preBötC influences arousal, cognitive function, and emotion (Del Negro et al., 2018). In turn, brain states (such as anxiety) lead to changes in respiration. However, even if these structures and communication between them are included in the respiratory modulation of the alpha CFPS(t) oscillations, at this stage, it remains unclear through which particular mechanism could respiration, assuming our hypothesis is valid, be influencing the observed CF phase shifts.

In addition to neurons in the PreBötzinger complex, respiration is actively influenced by groups of neurons in the pons and medulla, as well as central and peripheral chemoreceptors, mechanoreceptors, and metaboreceptors. Pontine neuronal groups incorporate the pneumotaxic and apneustic centers, which participate in control of the speed of breathing and perform fine-tuning of the rate of respiration (Webster and Karan, 2020). Dorsal medullar neurons are responsible for inhalation and airway defense, while the ventral neurons control exhalation (Bolser et al., 2015; Pagliardini et al., 2011). However, although mutually well-connected and functionally well-organized, the main external connection of these medullar groups with cortical structures is achieved through the PreBötzinger complex and locus coeruleus (Del Negro et al., 2018; Webster and Karan, 2020). This fact allows us to further hypothesize that their role in the emergence of the detected CFPS(t) rhythm is indirect. In contrast, sensory receptors project to the breathing centers not the cortex itself (Berger and Averill, 1983; Del Negro et al., 2018), so their role could also be considered indirect. Future investigations will determine whether and to what extent any of the abovementioned structures participate in generating and/or sustaining the detected 0.2232 Hz rhythm. At this stage, it is worth noting that, despite appearing in the cortex, the new phase shift rhythm was detected during spontaneous breathing, implying that the breathing process was not influenced by an active and conscious subject. How and whether this rhythm will be modified during paced breathing remains an important question to be answered in our future experiments.

A major limitation in the present work is the fact that we did not record both respiratory and EEG signals from the same subjects simultaneously and consequently have not been able to prove that the detected CFPS(t) oscillation at ≈ 0.22 Hz was phase-locked to the respiration. Understandably, this is our first objective for future research. However, according to our view, the reasonable congruence of the recorded frequencies between the two subject groups participating in two different sets of experiments, as shown in Figures 9C, D, leaves little doubt that the detected CFPS(t) oscillation originates from the subjects’ respiratory activity. However, respecting scientific rigor, we abstained from naming it “respiratory rhythm” but rather used the term “respiratory frequency rhythm” in the title and throughout this work.

 4.3 Alpha activity and respiration in different states of vigilance

If we compare the topographic distribution of five channel pairs with maximal FFT (CFPS(t)) oscillation amplitudes at 0.2232 Hz in wake and drowsy states (Figures 10B, 12B), one can see that in the wake state, they are predominantly longitudinal and laterally asymmetrical, positioned on the right side of the scalp. This particular distribution indicates that neural networks and populations positioned under and/or connecting anterior (frontal, F3, F4, and F8) and right posterior (P4, T6, and O2) scalp areas act in an organized way (underlying mechanisms unknown and remain yet to be elucidated) to participate in the respiratory frequency modulation of the detected alpha phase shifts. Considering the fact that neuronal excitability changes in synchrony with the instantaneous alpha phase (Romei et al., 2008; Sauseng et al., 2009; Mathewson et al., 2009; Berger et al., 2014), it follows that, for each channel pair, time delay of the excitation degree between these two large neuronal ensembles varies according to the respiratory frequency rhythm. This time delay, Δt, depends on the carrier frequency of the alpha activity, f ∝ (≈10 Hz), and the instantaneous value of CFPS(t), Δφ: Δ t = Δ φ 1 f ∝ 360 ° . Therefore, for each 10° of CFPS, Δt ≈ 2.78 m (compared to the CFPS waveform shown in the inlet of Figure 9C). The physiological significance of this rhythmically varying time delay is also not clear, but it obviously involves a more complex underlying mechanism related to, among others, the state of alertness and cognitive processes. However, since subjects did not perform any particular tasks, the topographic distribution of this rhythmicity in neuronal excitability may be more related to a state of readiness for incoming tasks, a kind of “phase shift resting state.” What appeared to be evident was that the large-scale respiratory frequency rhythmical integration of neuronal ensembles, involving anterior and posterior regions, seemed to exist in the wake state (Figures 10A, B). However, this picture was significantly altered in the drowsy state, where frontal regions appeared not to be rhythmically coupled with the posterior regions (Figures 12A, B). In terms of rhythmicity, anterior (frontal) and posterior regions became disentangled. The role of the frontal and prefrontal cortices across different states of alertness was already known (Muzur et al., 2002; Brown et al., 2012). Muzur et al. stated in their study that “functions of the prefrontal cortex are more relevant to the self-conscious awareness that is lost during sleep.” Similar results, pointing to the biggest amount of changes occurring in the frontal regions, were obtained in our previous analyses performed on the same set of signals during the wake to drowsy transition (Kalauzi et al., 2012b; 2015). Figures 10C, 12C present an additional result regarding the difference between the two states of vigilance: amplitudes of respiratory frequency CFPS(t) oscillations are greater in the wake state than those in the drowsy state, and this is valid for all channel pairs. Therefore, full wakefulness appeared to be necessary for this alpha phase shift oscillations to appear in their full capacity.

 4.4 Volitional control of breathing and alpha phase shifts

It was found that, besides the cortex, breathing rates can be influenced by stimulation of the hippocampus, amygdala, and insula (Frysinger and Harper, 1990; Harper et al., 1998; Herrero et al., 2018), indicating a complex network of different brain structures during the volitional control of breathing. In contrast, scalp EEG (Fumoto et al., 2004), as well as TMS (Locher et al., 2006), and neuroimaging techniques (McKay et al., 2003) revealed the role of premotor, motor, and supplementary motor cortices in the process. All results presented in our work were obtained in the regime of spontaneous breathing. Within this experimental framework and at this moment, it is not possible to predict how these results, specifically profiles (oscillatory content) of CFPS(t) spectra and topographic distribution of their most prominent amplitude peaks, would appear within a paced breathing paradigm. Potential appearance of peaks at the frequencies of paced breathing would certainly strengthen the hypothesis that it is linked to the breathing process and point to its possible role in the volitional control. Nevertheless, in order to elucidate the effect of volitional control, future experiments using topographic analyses of CFPS(t) respiratory frequency amplitudes under paced breathing should take care to place scalp electrodes over the cortical structures mentioned above.

Regarding the results of statistical tests presented in Figure 13, we found that, although amplitudes of alpha activity are reduced in the drowsy state compared to the wake state (Bojić et al., 2010), CFPS(t) oscillations did not disappear but still persisted at the frequency of 0.2232 Hz (although also reduced; Figures 10C, 11D, 12C). However, let us point out that we obtained an unexpectedly high number of channel pairs (20/91), which exhibited respiratory frequency amplitudes significantly larger in wake than in the drowsy state. To have 22% of such pairs, with only 10 subjects analyzed, seems encouraging and raises expectations that future experiments, if performed on a wider group of subjects, might produce more information and offer new insights into this newly discovered oscillation. For example, since our results point to a ≈0.2-Hz rhythmical variability of time delays between occurrences of extremal (minimal or maximal) neuronal excitability between different scalp regions, with simultaneously recorded EEG and respiratory signals, we would be able to determine whether these delays are (or not) locked to extremal values of the neuronal excitability caused by respiration, allowing us a clearer and more detailed picture of the interaction between respiration and alpha waves. In addition, a high percentage of channels with significantly different FFT respiratory frequency amplitudes between the two states of vigilance points to an important role of wakefulness and conscious awareness in the dynamical organization of respiratory modulation of alpha brain activity.

 4.5 Network approach: how to connect the brain, respiration, and cardiac activity

Cardiac rhythm and its causal relationship with EEG were studied in different sleep stages by different methodological approaches (e.g., by Granger causality; Abdalbari et al., 2022). In the course of this study, we searched for the presence of cardiac-like rhythm (around 1 Hz) in our FFT (CFPS(t)) spectra but could not detect it in any pair of electrodes. In contrast, in the case of simultaneous EEG and respiratory recordings, CFPS(t) time series may be included into multivariate predictability analyses such as those performed by Faes et al. (2016) in different sleep or pathological experimental paradigms (apnea-related conditions are particularly interesting states to explore).

Another set of situations where phase relations between respiration and CFPS(t) could be studied comprises different emotional psychological states since they involve both the brain and cardiovascular systems. It was shown that emotional response is initially caused by sympathovagal activity, in which ascending modulations precede efferent information transfer and correlate with the level of arousal (Candia-Rivera et al., 2022).

In course of the brain–body interaction studies, a network approach was employed involving multivariate correlation analysis (Pernice et al., 2021) and information dynamics (information storage and information transfer; Zanetti et al., 2019). More relevant to us was the work of Pernice et al., who calculated multivariate correlations within and between the human brain–body subnetworks in three experimental conditions representing different levels of mental stress: rest, sustained attention, and mental arithmetic. The quantities of interest to our present approach were the heart rate (η), respiration amplitude (ρ), and amplitudes of four EEG bands: δ, θ, α, and β. The authors found that internal subnetwork linkages were most common (mainly η–ρ; δ–θ; θ–α; and α–β). However, two inter-subnetwork correlations were also significant (η–β and η–δ). Therefore, when only amplitudes were being studied, neither alpha nor respiration was found to be mutually correlated. It would be an interesting alternative to probe whether alpha phases (phase shifts) are correlated with respiration phases in these states of heightened alertness as we found a considerable topographic redistribution of channels with maximal alpha FFT (CFPS(t)) peaks at 0.2232 Hz in the state of reduced alertness (drowsiness), compared to the nominal state (resting wakefulness).

Spontaneous breathing induces diverse autonomic rhythms that modulate cardiovascular control including arterial blood pressure and cerebral blood flow—heart rate and stroke volume— which, in turn, can influence EEG oscillations (Rau et al., 1993; Liu et al., 2019). Respiratory sinus arrhythmia, which has been intensively studied for decades (Hirsch and Bishop, 1981; Larsen et al., 2010), is an important part of this rhythmicity. All of the abovementioned rhythms initially originate from the breathing process itself, although the situation is more complex, as studies of causal relationships between cardiovascular and respiratory systems revealed a bidirectional influence (Tzeng et al., 2003; Cairo et al., 2023; Platiša et al., 2023; Porta et al., 2024). Therefore, having this network complexity in mind, if our hypothesis were to be proven, it would be a challenge to explore the role of RSA in the causal chain of these rhythmic processes. One possible approach might be to search for the link/correlation between the amplitude of RSA and FFT (CFPS (t)) peak height at ≈ 0.22 Hz in various experimental conditions and different physiological states.

Cardiorespiratory variables like the respiratory rate, heart rate, systolic arterial pressure, mean arterial pressure, and mean cerebral blood flow are main candidates for mediators by which the brain and respiration exert mutual influence. For instance, it has been shown that respiratory activity can act as a confounder or suppressor of the causal relationship between different cardiovascular and cerebrovascular variables. By applying the confounding/suppression test on propofol-based general anesthesia patients and controls, Porta et al. (2022) concluded that respiration behaved systematically as a confounder for cardiovascular and cerebrovascular controls. Causal relationships between the abovementioned variables can be quantitatively assessed by applying well-established methods both in time (Granger, 1980) and frequency (Geweke, 1982) domains. By studying postural syncope propensity, Porta et al. (2023) found that the postural syncope is favored by a loss of coordination between the baroreflex feedback and mechanical feedforward pathway in response to head-up tilt (HUT) in the LF band and by a weaker ability of cerebral autoregulation to limit mean cerebral blood velocity variability driven by mean arterial pressure changes in the respiratory rate during HUT. Having these results in mind, it cannot be a priori excluded that the respiratory impact on alpha phase shifts is mediated through pressure-to-flow at the cerebral level and baroreflex at the respiratory rate. Moreover, the relationship between the cerebral blood flow and alpha EEG activity (especially their phases; see Jann et al., 2009) must be modeled quantitatively in order to complete this cause–effect chain. Hopefully future analyses, specifically Geweke spectral causality applied to CFPS(t) and other abovementioned simultaneously recorded variables (primarily respiration), in spontaneous and paced breathing, will shed some new light on the causal relationships in the respiratory frequency range within and between respiratory, cardiovascular, and cerebrovascular systems. A carefully designed experimental setup, able to assess possible concomitant influences from respiration to EEG alpha activity (and vice versa), by mean arterial pressure to mean cerebral blood flow and to the heart period through baroreflex, at the frequencies of respiration and heart rate (respectively), is necessary.

 4.6 Potential future use of the detected rhythm

It is known that changes in cardiorespiratory coupling (CRC) are most expressed in stressful conditions. In some pathological states, such as type 2 diabetes mellitus, it was shown that they could be used as an early marker of cardiac autonomic dysfunction (Da Silva et al., 2023). The other set of conditions where CRC becomes increasingly strained involve physical exercise and sports training, where universal methodological standards have not yet been fully established (de Abreu et al., 2023; de Abreu et al., 2024). If our hypothesis is to be confirmed, we would have a newly detected EEG rhythm, probably more firmly correlated with respiration than existing EEG rhythms, in addition to conventional EEG–ECG–respiration recordings. In that case, present CRC could be extended to cardiorespiratory–brain coupling (CRBC), making it more sensitive to the changing conditions or pathological situations such as tumors, multiple sclerosis, or even degenerative diseases. We expect that in future experiments with CFPS(t)–ECG–respiratory signals, additional information could be acquired to enable us to reach those standards more easily.

 Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

 Ethics statement

The studies involving humans were approved by the Ethical Committee of the Faculty of Medicine, University of Belgrade (No. 2650/IV-24). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

 Author contributions

AK: conceptualization, data curation, formal analysis, methodology, software, validation, writing–original draft, and writing–review and editing. ZM: data curation, formal analysis, investigation, resources, supervision, writing–original draft, and writing–review and editing. ES: funding acquisition, project administration, supervision, visualization, writing–original draft, and writing–review and editing. TB: data curation, funding acquisition, investigation, resources, supervision, validation, visualization, writing–original draft, and writing–review and editing.

 Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This research was supported by Annex contract for institutional financing 451-03-9/2021-14/200017 of the Ministry of Science, Technological Development and Innovation of the Republic of Serbia and the project Proof of Concept 5537 supported by the Innovation Fund of the Republic of Serbia.

The authors express their gratitude to the Institute for Mental Health, Belgrade, and Laboratory for Biosignals, Institute for Biophysics, Faculty of Medicine, University of Belgrade, where part of this work was done, as well as to all kind women and men participating patiently in our experiments.

 Conflict of interest

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

 Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Abdalbari H. Durrani M. Pancholi S. Patel N. Nasuto S. J. Nicolaou N. (2022). Brain and brain-heart Granger causality during wakefulness and sleep. Front. Neurosci. 16, 927111. 10.3389/fnins.2022.927111 Bai X. Li J. Zhou L. Li X. (2009). Influence of the menstrual cycle on nonlinear properties of heart rate variability in young women. Am. J. Physiol. Heart Circ. Physiol. 297 (2), H765–H774. 10.1152/ajpheart.01283.2008 Berger A. J. Averill D. B. (1983). Projection of single pulmonary stretch receptors to solitary tract region. J. Neurophysiol. 49 (3), 819–830. 10.1152/jn.1983.49.3.819 Berger B. Minarik T. Liuzzi G. Hummel F. C. Sauseng P. (2014). EEG oscillatory phase‐dependent markers of corticospinal excitability in the resting brain. Biomed. Res. Int. 2014 (1), 936096. 10.1155/2014/936096 Bing-Canar H. Pizzuto J. Compton R. J. (2016). Mindfulness-of-breathing exercise modulates EEG alpha activity during cognitive performance. Psychophysiology 53, 1366–1376. 10.1111/psyp.12678 Bojić T. Vučković A. Kalauzi A. (2010). Modeling EEG fractal dimension changes in wake and drowsy states in humans – a preliminary study. J. Theor. Biol. 262 (2), 214–222. 10.1016/j.jtbi.2009.10.001 Bolser D. C. Pitts T. E. Davenport P. W. Morris K. F. (2015). Role of the dorsal medulla in the neurogenesis of airway protection. Pulm. Pharmacol. Ther. 35, 105–110. 10.1016/j.pupt.2015.10.012 Brændholt M. Kluger D. S. Varga S. Heck D. H. Gross J. Allen M. G. (2023). Breathing in waves: understanding respiratory-brain coupling as a gradient of predictive oscillations. Neurosci. Biobehav. Rev. 152, 105262. 10.1016/j.neubiorev.2023.105262 Brown R. E. Basheer R. McKenna J. T. Strecker R. E. McCarley R. W. (2012). Control of sleep and wakefulness. Physiol. Rev. 92, 1087–1187. 10.1152/physrev.00032.2011 Cairo B. Bari V. Gelpi F. De Maria B. Porta A. (2023). Assessing cardiorespiratory interactions via lagged joint symbolic dynamics during spontaneous and controlled breathing. Front. Netw. Physiol. 3, 1211848. 10.3389/fnetp.2023.1211848 Callaway E. Yeager C. L. (1960). Relationship between reaction time and electroencephalographic alpha phase. Science 132 (3441), 1765–1766. 10.1126/science.132.3441.1765 Candia-Rivera D. Catrambone V. Thayer J. F. Gentili C. Valenza G. (2022). Cardiac sympathetic-vagal activity initiates a functional brain–body response to emotional arousal. Proc. Nat. Acad. Sci. 119 (21), e2119599119. 10.1073/pnas.2119599119 Damoiseaux J. S. Rombouts S. A. Barkhof F. Scheltens P. Stam C. J. Smith S. M. (2006). Consistent resting-state networks across healthy subjects. Proc. Natl. Acad. Sci. U. S. A. 103 (37), 13848–13853. 10.1073/pnas.0601417103 Da Silva C. D. Catai A. M. de Abreu R. M. De Favari Signini É. Galdino G. A. M. Lorevice L. (2023). Cardiorespiratory coupling as an early marker of cardiac autonomic dysfunction in type 2 diabetes mellitus patients. Respir. Physiol. Neurobiol. 311, 104042. 10.1016/j.resp.2023.104042 Da Silva F. L. Vos J. E. Mooibroek J. Van Rotterdam A. (1980). Relative contributions of intracortical and thalamo-cortical processes in the generation of alpha rhythms, revealed by partial coherence analysis. Electroencephalogr. Clin. Neurophysiol. 50 (5–6), 449–456. 10.1016/0013-4694(80)90011-5 de Abreu R. M. Cairo B. Porta A. (2023). On the significance of estimating cardiorespiratory coupling strength in sports medicine. Front. Netw. Physiol. 2, 1114733. 10.3389/fnetp.2022.1114733 de Abreu R. M. Cairo B. Rehder-Santos P. da Silva C. D. De Favari Signini É. Milan-Mattos J. C. (2024). Cardiorespiratory coupling is associated with exercise capacity in athletes: a cross-sectional study. Respir. Physiol. Neurobiol. 320, 104198. 10.1016/j.resp.2023.104198 Del Negro C. A. Funk G. D. Feldman J. L. (2018). Breathing matters. Nat. Rev. Neurosci. 19 (6), 351–367. 10.1038/s41583-018-0003-6 Eke A. Hermann P. Bassingthwaighte J. B. Raymond G. M. Percival D. B. Cannon M. (2000). Physiological time series: distinguishing fractal noises from motions. Pflugers Arch. – Eur. J. Physiol. 439, 403–415. 10.1007/s004249900135 Ergenoglu T. Demiralp T. Bayraktaroglu Z. Ergen M. Beydagi H. Uresin Y. (2004). Alpha rhythm of the EEG modulates visual detection performance in humans. Cogn. Brain Res. 20, 376–383. 10.1016/j.cogbrainres.2004.03.009 Esteller R. Vachtsevanos G. Echauz J. Litt B. (2001). A comparison of waveform fractal dimension algorithms. IEEE Trans. Circuits Syst. I, Fundam. Theory Appl. 48 (2), 177–183. 10.1109/81.904882 Faes L. Marinazzo D. Stramaglia S. Jurysta F. Porta A. Giandomenico N. (2016). Predictability decomposition detects the impairment of brain–heart dynamical networks during sleep disorders and their recovery with treatment. Phil. Trans. R. Soc. A 374, 20150177. 10.1098/rsta.2015.0177 Fleming S. Thompson M. Stevens R. Heneghan C. Plüddemann A. Maconochie I. (2011). Normal ranges of heart rate and respiratory rate in children from birth to 18 years of age: a systematic review of observational studies. Lancet 377, 1011–1018. 10.1016/S0140-6736(10)62226-X Frysinger R. C. Harper R. M. (1990). Cardiac and respiratory correlations with unit discharge in epileptic human temporal lobe. Epilepsia 31 (2), 162–171. 10.1111/j.1528-1167.1990.tb06301.x Fumoto M. Sato-Suzuki I. Seki Y. Mohri Y. Arita H. (2004). Appearance of high-frequency alpha band with disappearance of low-frequency alpha band in EEG is produced during voluntary abdominal breathing in an eyes-closed condition. Neurosci. Res. 50, 307–317. 10.1016/j.neures.2004.08.005 Geweke J. (1982). Measurement of linear dependence and feedback between multiple time series. J. Am. Stat. Assoc. 77, 304–313. 10.2307/2287238 Granger C. W. J. (1980). Testing for causality. J. Econ. Dyn. Control 2, 329–352. 10.1016/0165-1889(80)90069-X Halgren M. Ulbert I. Bastuji H. Fabó D. Erőss L. Rey M. (2019). The generation and propagation of the human alpha rhythm. Proc. Natl. Acad. Sci. U. S. A. 116 (47), 23772–23782. 10.1073/pnas.1913092116 Harper R. M. Poe G. R. Rector D. M. Kristensen M. P. (1998). Relationships between hippocampal activity and breathing patterns. Neurosci. Biobehav. Rev. 22 (2), 233–236. 10.1016/S0149-7634(97)00010-9 Harsh J. Voss U. Hull J. Schrepfer S. Badia P. (1994). ERP and behavioral changes during the wake/sleep transition. Psychophysiology 31, 244–252. 10.1111/j.1469-8986.1994.tb02213.x Heck D. H. McAfee S. S. Liu Y. Babajani-Feremi A. Rezaie R. Freeman W. J. (2017). Breathing as a fundamental rhythm of brain function. Front. Neural Circuits 10, 115. 10.3389/fncir.2016.00115 Herrero J. L. Khuvis S. D. Yeagle E. Cerf M. Mehta A. D. (2018). Breathing above the brain stem: volitional control and attentional modulation in humans. J. Neurophysiol. 119, 145–159. 10.1152/jn.00551.2017 Hirsch J. A. Bishop B. (1981). Respiratory sinus arrhythmia in humans: how breathing pattern modulates heart rate. Am. J. Physiol. 241 (4), H620–H629. 10.1152/ajpheart.1981.241.4.H620 Hori T. Hayashi M. Morikawa T. (1990). “Topography and coherence analysis of the hypnagogic EEG,” in Sleep’90. Editor Horn J. (Bochum: Pontenagel Press), 10–12. Hsu S.-M. Tseng C.-H. Hsieh C.-H. Hsieh C.-W. (2020). Slow-paced inspiration regularizes alpha phase dynamics in the human brain. J. Neurophysiol. 123, 289–299. 10.1152/jn.00624.2019 Ito J. Roy S. Liu Y. Cao Y. Fletcher M. Boughter J. D. (2014). Whisker barrel cortex delta oscillations and gamma power in the awake mouse are linked to respiration. Nat. Commun. 5, 3572. 10.1038/ncomms4572 Jann K. Dierks T. Boesch C. Kottlow M. Strik W. Koenig T. (2009). BOLD correlates of EEG alpha phase-locking and the fMRI default mode network. NeuroImage 45, 903–916. 10.1016/j.neuroimage.2009.01.001 Javorka M. El-Hamad F. Czippelova B. Turianikova Z. Krohova J. Lazarova Z. (2018). Role of respiration in the cardiovascular response to orthostatic and mental stress. Am. J. Physiol. Regul. Integr. Comp. Physiol. 314, R761–R769. 10.1152/ajpregu.00430.2017 Jung F. Yanovsky Y. Brankačk J. Tort A. B. L. Draguhn A. (2023). Respiratory entrainment of units in the mouse parietal cortex depends on vigilance state. Pflugers Arch. – Eur. J. Physiol. 475, 65–76. 10.1007/s00424-022-02727-2 Jung J.-Y. Kang C.-K. (2021). Investigation on the effect of oral breathing on cognitive activity using functional brain imaging. Healthc 9, 645. 10.3390/healthcare9060645 Kalauzi A. Bojić T. Vuckovic A. (2012a). Modeling the relationship between Higuchi's fractal dimension and Fourier spectra of physiological signals. Med. Biol. Eng. Comput. 50 (7), 689–699. 10.1007/s11517-012-0913-9 Kalauzi A. Janković B. Ćulić M. Šaponjić J. Rajšić N. Šuljagić S. (1998). Phase demodulation of EEG signals. XLII ETRAN Conf. Proc. IV, 204–207. Kalauzi A. Kesic S. Saponjic J. (2009). Cortico-pontine theta synchronization phase shift following monoaminergic lesion in rat. J. Physiol. Pharmacol. 60 (4), 79–84. Kalauzi A. Matić Z. Bojić T. Platiša M. M. (2023a). Structure of Poincaré plots revealed by their graph analysis and low pass filtering of the RRI time series. Biomed. Signal Process. Control 80 (3), 104352. 10.1016/j.bspc.2022.104352 Kalauzi A. Matić Z. Platiša M. M. Bojić T. (2023b). Two operational modes of cardio-respiratory coupling revealed by pulse-respiration quotient. Bioeng. (Basel) 10 (2), 180. 10.3390/bioengineering10020180 Kalauzi A. Spasić S. Ćulić M. Grbić G. Martać Lj. (2005). Consecutive differences as a method of signal fractal analysis. Fractals 13 (4), 283–292. 10.1142/S0218348X05002933 Kalauzi A. Vuckovic A. Bojić T. (2012b). EEG alpha phase shifts during transition from wakefulness to drowsiness. Int. J. Psychophysiol. 86 (3), 195–205. 10.1016/j.ijpsycho.2012.04.012 Kalauzi A. Vuckovic A. Bojić T. (2015). Topographic distribution of EEG alpha attractor correlation dimension values in wake and drowsy states in humans. Int. J. Psychophysiol. 95 (3), 278–291. 10.1016/j.ijpsycho.2014.11.008 Kalauzi A. Vuckovic A. Bojić T. (2018). New complexity measures reveal that topographic loops of human alpha phase potentials are more complex in drowsy than in wake. Med. Biol. Eng. Comput. 56 (6), 967–978. 10.1007/s11517-017-1746-3 Klimesch W. (1999). EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Res. Rev. 29, 169–195. 10.1016/S0165-0173(98)00056-3 Klimesch W. Sauseng P. Hanslmayr S. (2007). EEG alpha oscillations: the inhibition–timing hypothesis. Brain Res. Rev. 53 (1), 63–88. 10.1016/j.brainresrev.2006.06.003 Klimesch W. Schack B. Sauseng P. (2005). The functional significance of theta and upper alpha oscillations for working memory: a review. Exp. Psychol. 52, 99–108. 10.1027/1618-3169.52.2.99 Kluger D. S. Balestrieri E. Busch N. A. Gross J. (2021). Respiration aligns perception with neural excitability. Elife 10, e70907. 10.7554/eLife.70907 Kluger D. S. Gross J. (2020). Depth and phase of respiration modulate cortico-muscular communication. NeuroImage 222, 117272. 10.1016/j.neuroimage.2020.117272 Kluger D. S. Gross J. (2021). Respiration modulates oscillatory neural network activity at rest. PLoS Biol. 19 (11), e3001457. 10.1371/journal.pbio.3001457 Larsen P. D. Tzeng Y. C. Sin P. Y. Galletly D. C. (2010). Respiratory sinus arrhythmia in conscious humans during spontaneous respiration. Respir. Physiol. Neurobiol. 174 (1-2), 111–118. 10.1016/j.resp.2010.04.021 Liu X. Pu Y. Wu D. Zhang Z. Hu X. Liu L. (2019). Cross-frequency coupling between cerebral blood flow velocity and EEG in ischemic stroke patients with large vessel occlusion. Front. Neurol. 10, 194. 10.3389/fneur.2019.00194 Locher C. Raux M. Fiamma M. N. Morélot-Panzini C. Zelter M. Derenne J. P. (2006). Inspiratory resistances facilitate the diaphragm response to transcranial stimulation in humans. BMC Physiol. 6, 7. 10.1186/1472-6793-6-7 Mandelbrot B. B. van Ness J. W. (1968). Fractional brownian motions, fractional noises and applications. SIAM Rev. 10, 422–437. 10.1137/1010093 Mathewson K. E. Gratton G. Fabiani M. Beck D. M. Ro T. (2009). To see or not to see: prestimulus α phase predicts visual awareness. J. Neurosci. 29 (9), 2725–2732. 10.1523/JNEUROSCI.3963-08.2009 Matić Z. Kalauzi A. Moser M. Platiša M. M. Lazarević M. Bojić T. (2022). Pulse respiration quotient as a measure sensitive to changes in dynamic behavior of cardiorespiratory coupling such as body posture and breathing regime. Front. Physiol. 13, 946613. 10.3389/fphys.2022.946613 Matić Z. Platiša M. M. Kalauzi A. Bojić T. (2020). Slow 0.1 Hz breathing and body posture induced perturbations of RRI and respiratory signal complexity and cardiorespiratory coupling. Front. Physiol. 11, 24. 10.3389/fphys.2020.00024 McKay L. C. Evans K. C. Frackowiak R. S. Corfield D. R. (2003). Neural correlates of voluntary breathing in humans. J. Appl. Physiol. 95 (3), 1170–1178. 10.1152/japplphysiol.00641.2002 Muzur A. Pace-Schott E. F. Hobson J. A. (2002). The prefrontal cortex in sleep. Trends Cogn. Sci. 6 (11), 475–481. 10.1016/s1364-6613(02)01992-7 Naatanen R. Picton T. (1987). The N1 wave of the human electric and magnetic response to sound: a review and an analysis of the component structure. Psychophysiology 24 (4), 375–425. 10.1111/j.1469-8986.1987.tb00311.x Nakamura N. H. Oku Y. Fukunaga M. (2023). “Brain–breath” interactions: respiration-timing–dependent impact on functional brain networks and beyond. Rev. Neurosci. 35 (2), 165–182. 10.1515/revneuro-2023-0062 Nielsen T. Stenstrom P. Takeuchi T. Saucier S. Lara-Carrasco J. Solomonova E. (2005). Partial REM-sleep deprivation increases the dream-like quality of mentation from REM sleep and sleep onset. Sleep 28 (9), 1083–1089. 10.1093/sleep/28.9.1083 Nunez P. L. Wingeier B. M. Silberstein R. B. (2001). Spatial‐temporal structures of human alpha rhythms: theory, microcurrent sources, multiscale measurements, and global binding of local networks. Hum. Brain Mapp. 13 (3), 125–164. 10.1002/hbm.1030 Ogilvie R. D. Simons I. A. Kuderian R. H. MacDonald T. Rustenburg J. (1991). Behavioral, event-related potential, and EEG/FFT changes at sleep onset. Psychophysiology 28 (1), 54–64. 10.1111/j.1469-8986.1991.tb03386.x Ogilvie R. D. Wilkinson R. T. (1984). The detection of sleep onset: behavioral and physiological convergence. Psychophysiology 21, 510–520. 10.1111/j.1469-8986.1984.tb00234.x Pagliardini S. Janczewski W. A. Tan W. Dickson C. T. Deisseroth K. Feldman J. L. (2011). Active expiration induced by excitation of ventral medulla in adult anesthetized rats. J. Neurosci. 31 (8), 2895–2905. 10.1523/JNEUROSCI.5338-10.2011 Perl O. Ravia A. Rubinson M. Eisen A. Soroka T. Mor N. (2019). Human non-olfactory cognition phase-locked with inhalation. Nat. Hum. Behav. 3 (5), 501–512. 10.1038/s41562-019-0556-z Pernice R. Antonacci Y. Zanetti M. Busacca A. Marinazzo D. Faes L. (2021). Multivariate correlation measures reveal structure and strength of brain–body physiological networks at rest and during mental stress. Front. Neurosci. 14, 602584. 10.3389/fnins.2020.602584 Platiša M. M. Radovanović N. N. Pernice R. Barà C. Pavlović S. U. Faes L. (2023). Information-theoretic analysis of cardio-respiratory interactions in heart failure patients: effects of arrhythmias and cardiac resynchronization therapy. Entropy 25, 1072. 10.3390/e25071072 Porta A. de Abreu R. M. Bari V. Gelpi F. De Maria B. Catai A. M. (2024). On the validity of the state space correspondence strategy based on k-nearest neighbor cross-predictability in assessing directionality in stochastic systems: application to cardiorespiratory coupling estimation. Chaos 34, 053115. 10.1063/5.0192645 Porta A. Gelpi F. Bari V. Cairo B. De Maria B. Tonon D. (2022). Categorizing the role of respiration in cardiovascular and cerebrovascular variability interactions. IEEE Trans. Biomed. Eng. 69 (6), 2065–2076. 10.1109/TBME.2021.3135313 Porta A. Gelpi F. Bari V. Cairo B. De Maria B. Tonon D. (2023). Concomitant evaluation of cardiovascular and cerebrovascular controls via Geweke spectral causality to assess the propensity to postural syncope. Med. Biol. Eng. Comput. 61, 3141–3157. 10.1007/s11517-023-02885-0 Rau H. Pauli P. Brody S. Elbert T. Birbaumer N. (1993). Baroreceptor stimulation alters cortical activity. Psychophysiology 30 (3), 322–325. 10.1111/j.1469-8986.1993.tb03359.x Rihs T. A. Michel C. M. Thut G. (2007). Mechanisms of selective inhibition in visual spatial attention are indexed by alpha-band EEG synchronization. Eur. J. Neursci. 25 (2), 603–610. 10.1111/j.1460-9568.2007.05278.x Romei V. Brodbeck V. Michel C. Amedi A. Pascual-Leone A. Thut G. (2008). Spontaneous fluctuations in posterior alpha-band EEG activity reflect variability in excitability of human visual areas. Cereb. Cortex 18 (9), 2010–2018. 10.1093/cercor/bhm229 Sauseng P. Klimesch W. (2008). What does phase information of oscillatory brain activity tell us about cognitive processes? Neurosci. Biobehav. Rev. 32 (5), 1001–1013. 10.1016/j.neubiorev.2008.03.014 Sauseng P. Klimesch W. Gerloff C. Hummel F. C. (2009). Spontaneous locally restricted EEG alpha activity determines cortical excitability in the motor cortex. Neuropsychologia 47 (1), 284–288. 10.1016/j.neuropsychologia.2008.07.021 Stam C. J. van Cappellen van Walsum A. M. Micheloyannis S. (2002). Variability of EEG synchronization during a working memory task in healthy subjects. Int. J. Psychophysiol. 46 (1), 53–66. 10.1016/S0167-8760(02)00041-7 Stam C. J. van Dijk B. W. (2002). Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate datasets. Phys. D. 163 (3), 236–251. 10.1016/S0167-2789(01)00386-4 Tallon-Baudry C. Bertrand O. Delpuech C. Pernier J. (1996). Stimulus specificity of phase-locked and non-phase-locked 40 Hz visual responses in human. J. Neurosci. 16 (13), 4240–4249. 10.1523/JNEUROSCI.16-13-04240.1996 Tzeng Y. C. Larsen P. D. Galletly D. C. (2003). Cardioventilatory coupling in resting human subjects. Exp. Physiol. 88, 775–782. 10.1113/eph8802606 Väyrynen T. Helakari H. Korhonen V. Tuunanen J. Huotari N. Piispala J. (2023). Infra-slow fluctuations in cortical potentials and respiration drive fast cortical EEG rhythms in sleeping and waking states. Clin. Neurophysiol. 156, 207–219. 10.1016/j.clinph.2023.10.013 Viczko J. Sharma A. V. Pagliardini S. Wolansky T. Dickson C. T. (2014). Lack of respiratory coupling with neocortical and hippocampal slow oscillations. J. Neurosci. 34 (11), 3937–3946. 10.1523/JNEUROSCI.3581-13.2014 Vuckovic A. Radivojevic V. Chen A. C. Popovic D. (2002). Automatic recognition of alertness and drowsiness from EEG by an artificial neural network. Med. Eng. Phys. 24 (5), 349–360. 10.1016/S1350-4533(02)00030-9 Webster L. R. Karan S. (2020). The physiology and maintenance of respiration: a narrative review. Pain Ther. 9, 467–486. 10.1007/s40122-020-00203-2 Zanetti M. Faes L. Nollo G. De Cecco M. Pernice R. Maule L. (2019). Information dynamics of the brain, cardiovascular and respiratory network during different levels of mental stress. Entropy 21, 275. 10.3390/e21030275