In physiological networks an important problem is to distinguish direct from indirect links. Figure 1 depicts direct and indirect links in simple three-node networks. In all three cases, the node corresponding to the source signal z seems to influence the target node, i.e., signal x. While there is a direct link from z to x in subfigure (a), there is only an indirect link in subfigure (b), mediated by signal y. Subfigure (c) shows another form of an indirect link, where the link between z and x is purely due to the common influence of y on both signals. Note that a time lag, indicated by the operator L , is important only in case (c). If the time lag from y to z was longer than the lag from y to x (i.e., α > β), the indirect link would change its direction and point from x to z instead.

There are more complex and mixed cases e.g., direct and indirect links between the same two nodes, and coexisting links [Bartsch and Ivanov (2014)], but the ones shown in Figure 1 are the most basic setups, and in the following they will be used to test and compare G-causality and BPRSA. By studying results for modeled data, we show that G-causality is more appropriate to distinguish these setups within a certain range of detection limits.

A simplified definition of Granger causality is: “Variable z Granger-causes (G-causes) variable x if knowledge about z improves the forecast of x” [Granger (1969)]. This reflects our common understanding of cause and consequence: The cause must precede the consequence in time and if z has no effect on x, we do not call it causal. This idea is formalized in the framework of autoregressive (AR) processes.

Under fairly general conditions a random process can be described by an AR model of order p [see, e.g., Lütkepohl (2005)]. Consider two AR models of the time series x _t , one including and one excluding information on z _t , x t = ∑ i = 1 p ϕ i 1 x t − i + w t 1 ; STD w 1 = σ 1 , (1) x t = ∑ i = 1 p ϕ i 2 x t − i + ψ i 2 z t − i + w t 2 ; STD w 2 = σ 2 . (2) Equation 1 models the value of x at time t as a weighted average over its own past plus a white noise process w with zero mean and standard deviation (STD) σ ⁽¹⁾. The weighting factors ϕ _i can be obtained by minimizing the error term w ⁽¹⁾. The past of x carries information about its own future, but this information is necessarily incomplete due to the statistical nature of x. The better x can be described (forecast) from its own past, the lower the standard deviation σ ⁽¹⁾ of the residual w ⁽¹⁾.

Equation 2 additionally considers information about z. If this information helps to model (forecast) x, then the standard deviation σ ⁽²⁾ will be reduced compared to σ ⁽¹⁾. The G-value G _z→x is a measure for the improvement of the forecast of x by including z, and therefore a measure of causality, G z → x = ln σ 1 σ 2 . (3) G _z→x quantifies G-causality, however it remains unclear whether 1) an obtained G-value is significantly different from zero, and 2) the link is direct or indirect.

To resolve issue 1), one can test the null hypothesis that all coefficients ψ i ( 2 ) of Eq. 2 are practically zero, i. e., G-causality does not exist. This is done by assuming an F or chi-squared distribution, given that the estimators of the least square method coefficients are asymptotically normally distributed [Neusser (2011)], and estimating the probability of a non-zero mean value. Because no other variables (such as y) are taken into account, this is a pairwise analysis method, which can identify links from z to x in all three cases of Figure 1.

In order to differentiate between direct and indirect links, i.e., to resolve issue 2), a conditional analysis is necessary. Consider the two extended AR models, x t = ∑ i = 1 p ϕ i 3 x t − i + τ i 3 y t − i + w t 3   STD w 3 = σ 3 , (4) x t = ∑ i = 1 p ϕ i 4 x t − i + τ i 4 y t − i + ψ i 4 z t − i + w t 4 STD w 4 = σ 4 . (5) Equations 4, 5 are the same as Eqs. 1, 2, except that the past of y is added to both. Hence, σ ⁽⁴⁾ is lower than σ ⁽³⁾ if and only if z adds information that is not already provided by y. Therefore, the conditional analysis will not show y-conditional G-causality G z → x y = ln σ 3 σ 4 (6) for the indirect z → x links shown in Figures 1B,C, but only for the direct link in Figure 1A. Besides, cases (b) and (c) can be distinguished by a pairwise analysis of z and y (disregarding x). This idea can be extended to a set of arbitrarily many variables, but for the scope of this work three variables (time series x _t , y _t , and z _t ) are sufficient.

Stationarity is an important prerequisite for the AR framework, because the process’ characteristics (i.e., the coefficients ϕ _i , ψ _i , τ _i ) must not change over time. A stationary process is a process with a constant mean and a finite covariance function that is invariant to shifts in time. This requirement is problematic, because many physiological signals are inherently non-stationary [Ivanov et al. (1996); Goldberger et al. (2002)]. Here, we probe stationarity with the Augmented Dickey-Fuller (ADF) test, which is widely accepted [Paparoditis and Politis (2018)] and based on AR modeling, so that it operates in the same framework as G-causality analysis.

In summary: G-causality is based on the improvement of a forecast by including additional data from other signals, and it can only be applied to stationary data. It can be used to distinguish between the three setups shown in Figure 1, if pairwise analysis and conditional analysis are applied. With G-causality the existence of a causal link can be decided as yes/no question with a statistical test, and quantified with the G-value.

While the G-causality approach is sensitive to non-stationarities, the Phase Rectified Signal Averaging method [PRSA, Bauer et al. (2006b)] and its bivariate version [BPRSA, Schumann et al. (2008)] have been developed to study noisy, non-stationary signals. The methods are especially suitable for quasi-periodic time-series, where perturbations reset the signal phase at random times. The BPRSA approach can easily be extended to an arbitrary number of signals. The idea is to align windows of the target signal x that are in the same phase with respect to one or more trigger signals (y and z) and average over all these windows. The procedure is described in detail in Schumann et al. (2008), see also Bauer et al. (2006b,a, 2009); here we only provide a brief overview.

The easiest and standard way to define trigger events from signal z is to consider all positions, where the signal increases, i. e., z t ν > z t ν − 1 for a trigger event at time t _ν . We denote all m trigger points by t _ν , ν = 1, …, m. Any other criterion that returns a Boolean value (trigger event or no trigger event) is possible, including criteria that are based on multiple signals. Each trigger event at t _ν leads to an anchor point x t ν of the target signal. Then, windows of width 2L are chosen around each anchor point x t ν , x t ν − L , x t ν − L + 1 … , x t ν + L − 1 . (7) The resulting BPRSA function for a potential z → x link is the point-wise average of x _t in all of these m windows, BPRSA j = 1 m ∑ ν = 1 m x t ν + j , j = − L , − L + 1 , … , L − 1 . (8)

The choice of anchor points is supposed to guarantee that in each window the index j = 0 is at a similar phase of the physiological process. The average over all windows is an in-phase superposition and therefore insensitive to non-stationarities (that are slower than the time scale L) and artifacts.

If there is no relation between z and x, and m is large, the resulting BPRSA time series will be constant everywhere, however with statistical fluctuations. This also happens if the choice of trigger points is not appropriate to reveal a relation between the signals, because it does not reflect the underlying processes. Any significant deviation from a constant value for any BPRSA _j must be interpreted as a relation between z and x, but not necessarily in the sense of G-causality. A positive (negative) peak in the BPRSA time series indicates the positive (negative) influence of a trigger event in z on x, or—in case the peak is at a negative index j—from x on z. Thus, BPRSA yields temporal information on cause and effect, just like G-causality. However, we would like to note that—unlike G-causality analysis—BPRSA does not model the time series data in any way, but relies on an averaging procedure assuming that the selected trigger event criterion is suitable for the relation between the considered signals and that the Central Limit Theorem holds. Therefore, it is expected that longer data are needed for a reliable identification of causality relations with BPRSA.

Since no test for the statistical significance of such a relation has yet been proposed ¹ , we have studied and compared four tests. The null hypothesis is that there is no causality between z and x. If this holds, the trigger points are randomly distributed and each of the BPRSA _j values, Eq. 8, is an average over independent random numbers and thus normally distributed due to the central limit theorem.

The following statistical tests are considered and compared:

1) The one-sided Kolmogorov-Smirnov test [Massey (1951)] measures the difference between the probability density functions of the BPRSA and a normal distribution, providing a p-value for the null hypothesis.

In summary: The existence of causality in the sense of BPRSA can be tested by checking whether the BPRSA _j values are normally distributed. In addition, the peak height can provide quantitative information on the link strength. Compared to G-causality, this method is less sensitive to non-stationarities, and it is a model-free approach. However, BPRSA cannot distinguish direct from indirect links unless more evolved trigger criteria could be established for a conditional analysis.

Physiological time series were derived from polysomnography (PSG) measurements that were recorded in several European sleep laboratories between September 1997 and April 2000 as part of the EU-project SIESTA [Klosch et al. (2001)]. Before any analysis, we chose 36 young, healthy subjects with excellent signal quality (young control group—YC, aged 29 ± 6), 36 elderly, healthy subjects (elderly control group—EC, aged 51 ± 10) and 43 age-matched, elderly subjects with an apnea-hypopnea index (AHI) of at least 10 per hour (OSA group, aged 51 ± 9). AHI is the mean number of apnea and hypopnea events per hour when considering a full-night sleep. Genders are distributed approximately equally in the YC (17 male, 19 female) and EC group (18 male, 18 female), but the OSA group consists mostly of male participants (38 male, five female). We address this in the results section. For each subject, we derived:

• Instantaneous heart rate H as the inverse RR-interval, i.e., the time between two successive heart beats.

Subsequently, signals H and B were interpolated to 1 Hz resolution, for EEG α, averages over non-overlapping windows of one second length were taken. The resulting three time series were all sampled at 1 Hz, and were averaged to resolutions of 2, 5, 10, 15, 30, and 60 s (i.e., 0.5, 0.2, 0.1, 0.067, 0.033, and 0.017 Hz) for further analysis (see Figure 4).

For each block of 30 seconds, sleep stages were scored by a sleep technician on the basis of the PSG signal following the rules by Rechtschaffen and Kales [Hobson (1969)]: light sleep (stages 1 and 2—LS), deep sleep (stages 3 and 4—DS) and REM sleep. Each triple of time series (H, B and EEG α) was partitioned into patches of continuous sleep of the same sleep stage, typically several minutes, and then normalized to zero mean and unity variance.

Since stationarity is a crucial precondition for G-causality analysis, a stationarity test must precede further analysis. We based our algorithm on the widely accepted Augmented Dickey-Fuller Test [ADF test, Paparoditis and Politis (2018)]. On the one hand, only stationary patches can be used for G-causality analysis, on the other hand the percentage of usable patches should be maximized in order to obtain the best possible statistics. A trade-off can be achieved through variation of the model order. For each continuous patch of the same sleep stage, all three signals of the node triple (i.e., heart rate, respiration rate and EEG amplitude) are tested for stationarity on a model of order 5. If all of them are stationary, the G-value is calculated according to Eq. 6. If at least one is nonstationary, the same procedure is repeated with model order 4, then with model order 3. If the process can still not be modeled as stationary at model order 3, the patch is split in half and the same procedure is applied to both shorter patches, motivated by the fact that shorter time series are more likely to be sufficiently free of trends and variability in variance and autocorrelations. A stopping condition is set when the patches reach a length of less than six times the time series resolution (i.e., 90 s length is the lower limit for time series of 15 s resolution), because this is the minimum length required for an AR model of order five. If this condition is met, the data is considered nonstationary and discarded.

Of all the data that are potentially available, the percentage of stationary patches is shown in Figure 5. It does not make sense to analyze data for resolutions broader than 15 s as most of these patches do not contain stationary data. This is mainly because of two reasons: Firstly, there are too few sufficiently long sleep stage epochs for resolutions > 15 s, and secondly, windows with an equal amount of data points but lower temporal resolutions cover longer sleep episodes, which are less likely to be stationary (i.e., compare high resolution 1 s time series with low resolution 60 s data). In our datasets, the second reason is dominant.

Even though resolutions above 15 s are of interest for research on long-term correlations, they cannot be included because there is not enough stationary data. For each group the conditional G-values (Eq. 6) for each time series resolution and each sleep stage were averaged and weighted by the length of the patch that they were calculated from in order to account for the varying length of sleep stages. ²

Error bars were calculated using a bootstrap method [Efron and Tibshirani (1993)]: Out of all G-values that comprise a data point, a new set of G-values was randomly drawn. This set is as large as the original one, but the same G-value can be picked several times, while others are not included in a particular sample. Overall, 100 such samples were drawn and for each sample the (unweighted) mean was calculated. The standard deviation of these means is an estimate for the standard deviation of the actual mean, the standard error, and is plotted as error bar in Figure 6.

Each analyzed data patch consists of only tens to hundreds of data points, which is not enough to reliably detect weak links with the G-causality test (see Figure 2). Therefore, we compare our results with surrogate data that allow for an alternative way to validate G-causality: If the average G-value obtained for a particular group and pair of signals is different from the G-value of the surrogate data, G-causality can be assumed. For our surrogate analysis, the respiration rate, heart rate and EEG amplitude data were taken from three different subjects, so that no G-causality can be expected between any of the three signals [Toledo et al. (2002); Bartsch et al. (2007)]. An alternative method would be the adjusted-amplitude Fourier transform (AAFT) method [Theiler et al. (1992); Lavanga et al. (2020)].

The results for all groups and the surrogate data are shown in Figure 6. The figure shows two of the six combinations in each panel, grouped into pairs of signals for each sleep stage. The error bars represent the standard error calculated from the bootstrap procedure described above. Particularly for LS and REM, the difference in directionality increases for larger time windows (i.e., lower resolutions at 10 and 15 s), however, for 15 s the error bars are quite large in many cases. For this reason, time series with 10 s resolution were chosen to generate the physiological networks presented in Figure 7. Thus, the strength of the connection line between two nodes (“network link”) is proportional to the G-value at a resolution of 10 s. While Figure 7 presents less information than Figure 6, it is more helpful to find patterns and recognize important differences between the groups.

It is important to note that a time series resolution of 10 s does not mean that the relevant processes happen within 10 s. The model that is used to calculate the G-values comprises three to five past terms, which means that the processes occur in time windows of 30–50 s, but by default neglect causalities from processes with dynamics on time scales shorter than 10 s.

Subfigures (a), (b) and (c) in Figures 6, 7 show the results for the YC, EC and OSA group, respectively. The results for YC and EC are qualitatively very similar: The interaction between heart and respiration rate is symmetrical in light and deep sleep across all time scales and takes low to medium values (G < 0.04 at all times). This changes in REM sleep, where the causality from respiration to heart rate clearly dominates over the opposite direction, especially towards larger time scales. The G-value for the causality from respiration to heart rate at time scale 2 seconds is slightly increased (see Figures 6A,B only), which is an indicator of respiratory sinus arrhythmia (RSA), a well-known effect of modulation of the heart rate within the breathing cycle. During inspiration, the heart rate accelerates and it slows down during expiration [Angelone and Coulter (1964)]. This RSA peak disappears for REM sleep, which is in agreement with Bond et al. (1973), who describe “a total dissociation between respiration and rhythmic heart rate variability” during REM sleep. However, while RSA (which acts on relatively short time scales on the order of a breathing cycle) disappears in REM sleep, the causal relation from respiration to heart rate during REM sleep is shifted to longer time scales. The RSA peak is also absent in OSA subjects during non-REM sleep (see Figure 6C), possibly indicating reduced RSA for OSA subjects.

The coupling between EEG α amplitude and respiration rate remains, for all groups, constantly weak and symmetrical throughout all sleep stages and time series resolutions. The only deviation from this behavior is a slight increase in the causality from EEG α to respiration rate at larger time scales for EC and OSA subjects in REM sleep, which could be related to aging. Regarding the coupling between heart rate and EEG α, there is a clear dominance from EEG α to heart rate during light and REM sleep. This coupling almost completely vanishes during deep sleep, causing EEG α to heart rate coupling to become more symmetric.

The network plots 7(a) and (b) show the small differences between the YC and EC group. Furthermore, throughout all groups there is only little G-causality during deep sleep. This is in accordance with results of Bashan et al. (2012) and Bartsch et al. (2015) who show low network connectivity during deep sleep using a time delay stability approach to quantify interactions. At the same time, during deep sleep there is also a loss of long-term correlations in heart beats [Bunde et al. (2000); Penzel et al. (2003b)] as well as in respiratory inter-breath intervals [Schumann et al. (2010)]. Such long-term correlations, however, exist during light and REM sleep and are assumed to be due to influences from the sympathetic nervous system on cardiac and respiratory dynamics [Schmitt et al. (2009); Bashan et al. (2012)]. In contrast, deep sleep, which is considered the most restorative sleep stage [Dijk (2009)], is characterized by sympathetic withdrawal and greatly reduced influence of the autonomic nervous system on heart and respiratory dynamics. In our results, this is reflected by a more autonomous behavior of all three nodes for all groups. A careful comparison between the YC and EC group reveals that EEG α to heart rate coupling may slightly increase with age. For OSA patients this effect is even more pronounced and accompanied by an overall increase of coupling strength also between heart and breathing as well as EEG α to breathing, possibly indicating a decreased deep sleep quality because of sleep apnea.

In general, the OSA group follows similar sleep-stage patterns as the YC and EC groups. Distinct differences, however, can be seen in the overall strength of coupling: First and most strikingly, the G-causality from heart rate to respiration rate is much stronger during light and REM sleep than in the other two groups. This can be attributed to different relaxation speeds of heart and respiratory rate after apnea events. Following an episode of sleep apnea, heart and respiration rate are increased [Penzel et al. (2003a; 2016)], but the relaxation of the heart rate happens faster than the relaxation of the respiratory rate. This leads to a situation where changes in heart rate precede changes in respiratory rate and thus lead to a detection of increased G-causality, as indicated by higher G-values. Secondly, there is an increased G-causality from EEG α amplitude to breathing during REM sleep.

Figures 6D, 7D show results for the surrogate data. As expected, the G-value is very close to zero in all cases at all resolutions, with small deviations for the 15 s data point, which could be due to the large error bars at this resolution. These results can be used as a baseline showing that randomly assembled physiological data yield negligible G-values, and that any deviation from zero, as seen in all other sub-figures, are physiologically meaningful results.

Performing the same analysis without testing for stationarity yields slightly different results (see Supplementary Figures S4, S5). Values at short timescales are mostly unchanged, but it seems that the G-value is overestimated in some cases when ignoring the stationarity test. Additionally, the standard errors are larger, especially for longer time scales. Therefore, we conclude that only stationary data yield non-spurious network links when applying the G-causality method. The fact that the OSA group consists of mostly male subjects does not influence our results. We repeated the analyses of Figures 6, 7 excluding all female subjects from YC and EC, and obtained similar results (see Supplementary Figures S6, S7).